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ABSTRACT 


This  paper  deals  with  a  preliminary  evaluation  of  a  mean- 
squared  estimated  algorithm  for  the  extraction  of  a  signal  that 
is  embedded  in  superposed  stationary  noise  processes.  The  mo¬ 
tivation  for  this  study  is  the  processing  and  interpretation  of 
data  obtained  in  recent  experiments  to  measure  gradients  of 
magnetic  fields  induced  by  internal  waves  in  a  shallow  ocean. 

The  approach  used  is  simulation  of  the  signal  and  noise 
processes  through  construction  of  random  time  series  for  the 
signal  and  interfering  noise  processes  whose  spectra  can  readily 
be  adjusted  to  study  parametrically  a  wide  range  of  situations. 

The  principal  conclusions  are:  (1)  the  algorithm,  in  the 
form  commonly  applied  to  the  reconstruction  of  the  output  of  the 
main  data  channel,  is  valid  only  if  a  very  special  form  of  linear 
relationship  is  assumed  between  each  of  the  superposed  time  series 
in  the  main  channel  and  the  time  series  in  each  of  the  subsidiary 
channels;  an  assumption  that  can  lead  to  grossly  erroneous  esti¬ 
mates  of  the  signal  spectrum  and  (2)  severe  cumulative  bias  errors 
can  be  introduced  in  the  reconstructed  signal,  partial  coherences 
and  transfer  functions  through  the  algebraic  operations  necessary 
in  the  application  of  the  algorithm.  This  procedure  can  lead  to 
large  uncertainties  in  the  functional  form  of  the  reconstructed 
signal  spectra  and,  even  more  important,  to  false  indications  of 
the  state  of  partial  coherence. 
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SUMMARY  AND  CONCLUSIONS 


This  paper  deals  with  a  preliminary  evaluation  of  a  mean- 
squared  error  estimation  algorithm  for  the  extraction  of  a  sig¬ 
nal  that  is  embedded  in  superposed  stationary  noise  processes. 
The  signal,  together  with  several  additive  noise  processes, 
comprises  the  main  sensor  channel.  Statistical  estimates  of 
one  or  more  of  the  individual  noise  processes  as  well  as  of  the 
signal  are  also  available  from  several  independent  sensors 
(termed  subsidiary  channels),  which  data  is  employed  to  estimate 
the  relative  spectral  content  of  the  signal  and  noise  constitu¬ 
ents  in  the  main  channel.  The  principal  application  of  the  al¬ 
gorithm  has  been  in  the  validation  of  theoretical  response  char¬ 
acteristics  of  linear  sensors.  The  sensor  to  be  evaluated  is 
usually  taken  as  the  main  channel.  Subsidiary  channel  sensors 
are  assumed  to  possess  generally  known  response  characteristics 
and  to  provide  intrinsic  signal-to-instruiiient  noise  ratios  of 
sufficient  magnitude  to  yield  statistically  significant  infor¬ 
mation  on  the  principal  noise  constituents  in  the  main  channel. 

The  direct  motivation  for  this  study  is  the  processing 
and  interpretation  of  data  obtained  in  recent  experiments  to 
measure  gradients  of  magnetic  fields  induced  by  internal  wave 
activity  in  a  shallow  ocean.  The  experiment  has  been  carried 
out  off  the  coast  of  San  Diego  by  Physical  Dynamics,  Inc.  under 
ARPA  sponsorship.  The  purpose  of  the  experiment  was  to  validate 
the  theoretical  relationships  between  magnetic  field  gradients 
induced  by  internal  wave  motion  of  sea  water  relative  to  the 
geomagnetic  field  and  the  hydrodynamic  parameters  of  internal 
waves.  The  findings  are  expected  to  provide  data  for  the  design 
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of  future  at-sea  experiments,  the  objective  of  which  will  be 
the  verification  of  theoretically  predicted  ocean  wave  induced 
noise  constituents  in  superconducting  gradiometers  when  employed 
in  magnetic  anomaly  detection. 

Physical  Dynamics  has  used  the  aforementioned  algorithm  to 
reconstruct  the  measured  gradiometer  output  with  the  aid  of  data 
obtained  from  nine  subsidiary  sensors  (channels).  These  included 
current  flowmeters  and  thermistor  chains  yielding  data  on  the 
signal  (internal  waves)  as  well  as  sensors  monitoring  the  major 
main-channel  noise  constituents  (e.g.,  those  induced  by  random 
motion  of  the  measurement  platform).  Despite  the  fact  that  the 
algorithm  (in  essentially  the  same  form  as  employed  by  Physical 
Dynamics)  has  found  extensive  application  in  the  interpretation 
of  measurement  data,  a  detailed  evaluation  of  its  capabilities 
under  controlled  conditions  does  not  appear  to  have  been  under¬ 
taken. 

Because  the  interpretation  of  the  measurement  data  in  the 
Physical  Dynamics  experiment  as  well  as  in  other  current  and 
future  ARPA  nonacoustic  ASW  programs  appears  to  be  critically 
dependent  on  the  efficacy  of  this  algorithm,  IDA  was  directed 
to  undertake  a  limited  effort  to  investigate  the  algorithm's 
performance  characteristics  under  controlled  conditions,  using 
simulated  data.  In  this  preliminary  study,  no  attempt  has  been 
made  at  a  duplication  of  the  functional  forms  of  the  spectra 
observed  in  the  Physical  Dynamics  experiment.  Rather,  the 
approach  has  been  to  construct  random  time  series  for  the  sig¬ 
nal  and  interfering  noise  processes  whose  spectra  can  readily 
be  adjusted  to  encompass  a  sufficiently  wide  range  of  situations 
so  as  to  increase  the  likelihood  of  uncovering  the  major  intrin¬ 
sic  limitations  of  the  data  processing  procedure  per  se.  Thus, 
the  term  "gradiometer"  as  employed  in  the  sequel  is  to  be  taken 
merely  as  a  descriptive  designation  of  the  main  data  channel, 
and  is  not  to  be  construed  as  implying  any  relation  between  the 
simulated  spectra  and  those  that  might  be  observed  with  real 
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gradiometers .  Also,  instead  of  employing  the  complete  set  of 
nine  channels,  the  numerical  simulation  was  carried  out  only 
in  the  special  case  of  one  main  data  channel  and  two  subsidiary 
channels.  Nevertheless,  the  inherent  deficiencies  of  the  algo¬ 
rithm  as  enumerated  in  the  following  are  expected  to  apply  a 
fortiori  in  case  more  than  three  channels  are  employed.  The 
principal  conclusions  are  as  follows: 

•  The  form  of  the  algorithm  commonly  applied  to  reconstruct 
the  output  of  a  sensor  is  valid  only  if  a  very  special  type  of 
linear  relationship  subsists  between  each  of  the  hypothesized 
time  series  in  the  main  channel  and  the  time  series  in  each  of 
the  subsidiary  channels  (Eq.  15).  We  term  this  relation  func¬ 
tional  decoupling.  The  existence  of  functional  decoupling  would 
in  general  have  to  be  inferred  on  the  basis  of  a  priori  knowl¬ 
edge  since  its  direct  experimental  verification  is  not  always 
possible.  Moreover,  the  consequences  of  incorrectly  assuming 
functional  decoupling  can  lead  to  grossly  erroneous  estimates 

of  the  signal  spectrum  [Eq.  (65)]  versus  Eq.  (67)]. 

•  Because  algebraic  operations  currently  incorporated  in 
the  algorithm  to  determine  the  signal,  partial  coherences,  and 
transfer  functions  are  based  on  mathematical  results  derived 
for  ideal  spectra  (i.e.,  exact  Fourier  transforms  of  the  re¬ 
spective  correlation  functions),  cumulative  bias  errors  can 
result  when  the  same  operations  are  carried  out  with  biased 
estimators  of  these  spectra.  Although  the  latter  procedure 
conforms  to  standard  practice,  numerical  simulations  indicate 
that  it  can  lead  to  significant  uncertainties  in  the  recon¬ 
structed  signal  spectra  and,  what  is  even  more  important,  to 
false  indications  of  the  state  of  partial  coherence  between 
the  main  and  a  subsidiary  channel.  While  these  problems  may 
be  partially  alleviated  by  a  judicious  choice  of  data  windows, 
any  attempt  at  an  optimum  solution  would  of  necessity  involve 
a  reformulation  of  the  algorithm  directly  in  terms  of  spectral 
estimators . 
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•  Additional  observations  based  on  the  specific  parameters 
employed  in  the  simulation  study  may  be  found  under  ’’Concluding 
Remarks"  on  page  159- 
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1.  INTRODUCTION 


1.1  EXPERIMENTAL  BACKGROUND 

Seawater,  by  virtue  of  its  electrical  conductivity,  gener¬ 
ates  electrical  currents  through  its  motion  relative  to  the 
earth's  magnetic  field.  These  electric  currents  generate,  in 
turn,  magnetic  fields,  which  may  be  used  to  detect  internal 
wave  motions  below  the  surface.  Recent  measurements  of  magne¬ 
tic  field  gradients  by  Physical  Dynamics  (Ref.  1),  7m  above  the 
surface  of  l8m  deep  seawater,  appear  to  show  correlations  be¬ 
tween  magnetic  field  gradients  and  wave  motions  having  charac¬ 
teristic  frequencies  for  swell  (~100  mHz),  surfbeat  ( — 30  mHz) 
and  internal  waves  (~3  mHz).  The  measurements  were  performed 
at  an  oceanographic  research  tower  located  about  one  mile  off¬ 
shore  near  San  Diego,  California.  The  basic  instrumentation 
incorporated  nine  channels  to  (a)  measure  magnetic  gradients 
from  ocean  waves  with  emphasis  on  those  from  internal  waves 
and  (b)  monitor  magnetic  gradient  noise  from  conduction  and 
magnetization  currents  in  nearby  objects  as  well  as  instrument 
motions.  The  instrumentation  consisted  of  a  two-input  gradiom- 
eter  and  triaxial  fluxgate  magnetometer  above  the  sea  surface 
as  well  as  four  electromagnetic-type  biaxial  current  meters  at 
several  depths  below  the  sea  surface.  Other  supporting  instru¬ 
mentation  involved  two  biaxial  tilt  meters,  pressure  transducers, 
thermistors  and  other  sensors  giving  air  pressure,  temperature 
and  wind  speeds.  The  complete  data  records  included  fluctua¬ 
tions  of  the  ambient  magnetic  field  above  the  sea  surface,  sea¬ 
water  displacements  and  currents  below  the  sea  surface,  motions 
of  the  gradiometer,  as  well  as  the  ambient  pressure,  temperature 
and  wind  speeds. 
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The  above  measurements  provided  records  of  voltage  fluctu¬ 
ations  during  time  intervals  from  six  to  eight  hours  at  the  rate 
of  once  every  second.  To  preclude  aliasing,  the  data  were  first 
filtered  to  eliminate  frequencies  above  300  mHz.  The  data  rec¬ 
ords  were  then  averaged  every  12  seconds  to  obtain  each  point  of 
corresponding  time  series  of  2048  discrete  points.  Past  Fourier 
transform  (PFTs)  operations  on  these  time  series  gave  corres¬ 
ponding  auto-  and  cross-spectral  density  functions  having  1024 
frequencies  spaced  uniformly  over  a  frequency  range  from  0.41 
to  41.7  mHz.  Each  spectral  density  function  was  arithmetically 
averaged  over  32  neighboring  frequencies  so  as  to  reduce  their 
fluctuations  about  their  mean  values.  An  algorithm  for  a 
multiple-input  linear  system  was  then  utilized  to  obtain  esti¬ 
mates  of  (a)  the  transfer  functions  of  each  subsidiary  channel, 

(b)  the  spectrum  of  the  magnetic  field  gradients  or  desired  signal 
from  the  gradiometer  (Sg^^),  (c)  the  spectra  of  other  relevant 
parameters  including  those  of  the  air  pressure,  temperature  and 
wind  speeds,  and  (d)  the  relevant  partial  coherences.  A  sample 
estimate  of  the  mean  values  of  the  spectrum  for  the  magnetic 
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gradient,  Sg.^,  indicated  a  frequency  dependency  given  by 
Sgfg  =:  200f-2  within  the  frequency  range  0.4  <  f  <10  mHz, 
where  the  units  of  Sg^g^  are  given  in  (pT/m)2/Hz. 

1.2  theoretical  background 

The  central  problem  in  the  analysis  of  experimentally- 
obtained  data  records  (i.e.,  time  series)  is  the  extraction  of 
one  or  more  desirable  signals  (e.g.,  magnetic  gradients  genera¬ 
ted  by  internal  waves)  from  the  invariably  present  "noise"  back¬ 
ground  given  by  the  superposition  of  other  time  series  (e.g., 
conduction  and  magnetization  currents  other  than  those  from  in¬ 
ternal  waves  and  Instrument  motions).  In  many  instances,  the 
background  "noise"  is  decomposable,  at  least  in  part,  into  a  set 
of  distinct  time  functions  engendered  by  processes  with  known 
causative  mechanisms  susceptible  to  independent  subsidiary 
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measurements.  It  then  becomes  possible,  at  least  in  principle, 
to  employ  the  data  from  these  subsidiary  measurements  in  con¬ 
junction  with  a  suitable  estimation  algorithm  to  effect,  essen¬ 
tially,  an  enhancement  of  the  desired  signal-to-background  noise 
ratio  through  a  partial  subtraction  of  the  "noise"  time  series. 

If,  in  addition,  one  of  the  subsidiary  measurements  yields  data 
linearly  related  to  the  signal  itself,  an  algorithm  based  on  a 
suitable  linear  combination  of  all  the  subsidiary  data  channels 
can  then  be  used  to  provide  a  partial  reconstruction  of  the  main 
channel  data  record.  In  this  latter  case,  instead  of  noise  sub¬ 
traction,  one  considers  the  modeling  of  the  physical  processes 
that  give  rise  to  the  superposition  of  the  time  series  in  the 
main  channel.  In  contrast  with  the  signal-to-noise  enhancement 
of  the  main  channel  through  noise  subtraction,  the  reconstruc¬ 
tion  of  the  main-channel  data  record  constitutes  mostly  a  veri¬ 
fication  procedure  of  the  main-channel  data.  However,  such 
verifications  are  usually  only  partial,  because  not  all  the 
processes  contributing  to  the  time  series  in  the  main  channel 
are  susceptible  to  Independent  observations.  (For  example,  no 
independent  observations  of  intrinsic  instrument  noise  are  pos¬ 
sible  .  ) 

A  common  measure  of  the  fidelity  of  reconstruction  of  the 
several  time  series  in  the  main  channel  is  the  mean-squared 
error  between  the  difference  of  the  main  channel  output  and  a 
linear  superposition  of  all  the  outputs  of  the  subsidiary  chan¬ 
nels.  If  the  signal  and  "noise"  time  series  can  reasonably  be 
modeled  as  stationary  stochastic  processes,  the  algorithm  in¬ 
corporating  the  analytical  solution  for  the  minimization  of  the 
mean-squared  error  is  given  by  a  multi-channel  Wiener  Filter. 

The  actual  implementation  of  the  filter  can  be  of  either  the 
open  loop  or  adaptive  type.  Although  the  latter  has  several 
advantages,  only  the  former  Is  of  interest  herein.  Even  though 
this  type  of  algorithm  has  been  in  existence  for  over  thirty 
years,  it  is  somewhat  surprising  that  its  formulation  has  remained 
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limited  to  predicting  the  fidelity  of  reconstruction  under  only 
statistical  ensemble  averaging  or,  equivalently,  averaging  over 
very  large  (infinite)  time  records.  The  algorithm  yields,  for 
example,  the  transfer  function  of  each  channel  based  on  ideal 
auto-  and  cross-spectral  density  functions,  which  are  assumed 
to  be  determinable  from  ensemble  averaging  of  the  corresponding 
time  series. 

In  actual  practice,  the  algebraic  operations  dictated  by 
the  analytical  solution  under  assumptions  of  perfect  ensemble 
averaging  are  applied  to  statistical  estimates  of  the  self  and 
cross  spectra  obtained  via  FFT  operations  and  frequency  averag¬ 
ing.  Since  these  statistical  estimates  exhibit  fluctuations  as 
well  as  bias  errors,  the  final  quantities  of  interest  (e.g., 
reconstructed  signal,  partial  coherence)  can  deviate  substan¬ 
tially  from  the  "true"  statistical  averages. 

As  usual,  the  ideal  auto-spectral  density  function,  S-^Cf), 
of  the  time  series  x.^(t)  is  defined  by 

/°°  -12TrfT 

R11(t)e  dx  ,  (1) 

—  CO 

where  f  is  the  frequency  and  R-^(x)  denotes  the  auto-correlation 
function 

Rfl (t )  =  <x1(t+x)  x  (t)>  } 

wherein  the  notation  <  >  denotes  the  ensemble  or  statistical 
average  of  the  stationary  time  series  x1(t).  Similarly,  the 
ideal  cross-spectral  density  function  for  two  given  time  series 
is  given  by 

Vf>  'Ll  Riy(T,e'12”fTdT  (3) 

where  denotes  the  cross-correlation  function 
Rly(T)  =  <x1(t+x)  y(t)>  . 
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The  application  of  the  mean-squared  estimation  algorithm  utiliz¬ 
ing  the  above  ideal  spectral  density  functions  would  yield  the 
corresponding  ideal  transfer  function  for  each  channel  as  well 
as  the  ideal  mean  values  of  the  spectra  for  the  desired  signal. 
Since  under  the  usual  experimental  constraints  direct  use  of 
ensemble  averaging  is  impractical,  it  becomes  necessary  to  em¬ 
ploy  statistical  estimates  of  the  correlation  function  based 
on  time  averaging.  Since  the  advent  of  the  PFT  technique,  it 
has  become  common  practice  to  estimate  the  spectrum  directly 
in  the  frequency  domain  instead  of  first  computing  the  corre¬ 
sponding  time-averaged  correlation  functions.  Using  the  avail¬ 
able  single  and  finite  length  time  series  from  a  given  channel, 
the  estimated  auto-spectral  density  function,  for  example,  is 
obtained  by  computing  first  the  unaveraged  periodogram  of  the 
data  record  in  the  time  domain  0  s:  t  s  T  from 

Su(f,T)  s  f  \f0  x1(t)e‘i2lTftdt  (5) 

or 

Su(f,T)  -  X^f,T)  =  i  X^f.T)  2  ,  (6) 


where  X-^fjT)  is  the  finite  Fourier  transform  of  x^t),  i.e., 

T 

Xl(f,T)  =  /  xi(t)e'i2irftdt  (7 


ft 

and  X^CfjT)  is  the  complex  conjugate  of  X-^(f , t) .  Similarly,  the 
cross-spectral  density  function  for  a  pair  of  associated  sample 
records  x-^Ct)  and  y(t)  in  the  time  domain  0  <  t  <  T  is  given  by 


Siy(f.T) 


^  X*(f,T)  Y(f,T) 


(8) 


where 


y(t)e 


— 1  2tt  ft 


dt  . 


(9) 
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A  A 

The  unaveraged  periodogram  S^CfjT)  and  Sly(f,T)  are  determined 
via  PFTs  of  x.^(t)  and  y(t)  to  obtain  X^CfjT)  and  Y(f,T).  They 
are  usually  characterized  by  a  significant  variance,  i.e.,  fluc¬ 
tuations  about  their  mean  values.  If  the  observation  period 
T  were  sufficiently  large,  the  statistical  (ensemble)  mean  of 

A 

S,,(f,  T)  would  tend  to  S  , (f )  as  given  by  Eqs.  (1)  and  (2), 

±1.  -Li.  ^ 

i.e.,  to  the  spectrum  of  the  data  record  x^(fc).  However, 

(f,  T)  is  not  a  useful  estimator  because  its  variance  does  not 
diminish  regardless  of  how  large  T  is  chosen.  In  statistical 

A 

terminology,  the  estimator  Si;L(f,T)  as  given  by  Eq.  (5)  is  said 

to  be  inconsistent.  In  practice,  this  problem  is  circumvented 

by  introducing  a  second  step  to  smooth  out  the  fluctuations  of 
£ 

S-^CfjT)  by  using  average  values  at  each  frequency  over  several 
adjacent  frequencies.  For  example,  if  26  denotes  the  frequency 
interval  over  which  S11(f,T)  is  averaged  at  each  frequency, 
the  averaging  operation  may  be  represented  by* 

1  rf+ 6  * 

s11(f,6,T)  -  26  /  Su(n,T)dn  •  (10) 

J  f-6 

In  Its  discrete  form,  Eq.  (10)  involves  an  arithmetical  averaging 
over  M  frequencies  at  either  side  of  a  particular  frequency, 
i.e., 

f  +  6 

Sn(f,M,T)»  2MTI  n^6  Su(f,T,n)  (11) 

and  6  =  2irM/T. 

A 

The  estimator  of  the  spectrum,  (f,M,T)  is  usually  re¬ 

ferred  to  as  a  smoothed  periodogram.  A  desirable  consequence 

/\ 

of  frequency  averaging  is  that  the  variance  of  S ^  (f ,M,T)  can, 

Thedata  smoothing  and  the  Fourier  transforms  are,  of  course, 
performed  only  at  discrete  values  of  the  record.  This  aspect 
is  not  important  to  the  present  discussion.  For  convenience, 
we  assume  a  continuous  data  record  x1(t). 
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in  principle,  be  made  as  small  as  desired  by  choosing  6T  suffi- 

A 

ciently  large.  In  statistical  terminology,  the  estimator 
(f,M,T)  is  said  to  be  consistent.  Since  the  data  record  is 
usually  of  fixed  length,  the  reduction  in  variance  is  effected 
by  increasing  the  frequency  averaging  interval  26.  Unfortunately, 

A 

as  6T  is  increased,  the  statistical  mean  of  S^(f,M,T)  does  not 
tend  to  the  true  spectrum  of  x.,(t)  but  to  some  distorted  version 
of  the  spectrum.  The  estimator  S^(f,M,T)  is  then  said  to  be 
biased.  The  net  result  of  increasing  6  is  to  produce  estimators 
that  are  consistent  but  biased,  i.e.,  estimators  with  reduced 
fluctuations  about  mean  values  that  are  different  from  the  ideal 
or  true  spectrum,  the  difference  being  made  progressively  greater 
as  6  increases.  It  is  this  tradeoff  between  variance  and  bias 
with  increasing  6  that  represents  an  unresolved  thorny  problem 
in  the  interpretation  of  periodogram  averaged  spectra. 

The  problem  of  the  proper  interpretation  of  the  spectral 
features  is  magnified  substantially  when  use  is  made  of  smoothed 
periodograms  in  the  several  algebraic  operations  (e.g.,  matrix 
inversions)  required  by  the  algorithm  to  derive  the  quantities 
of  interest  such  as  transfer  functions  and  partial  coherences. 
Since  the  analytical  solution  for  the  bias  error  and  variance 
estimates  for  even  one  channel  is  extremely  difficult,  the 
most  appropriate  methodology  for  studying  these  effects  in  a 
multichannel  situation  is  numerical  simulation.  For  example, 
a  particular  stationary  process  with  a  known  spectrum  (signal) 
together  with  a  number  of  "noise"  spectra  can  readily  be  created 
by  appropriate  digital  filtering  of  random  number  sequences. 

The  spectral  reconstruction  algorithm  can  then  be  applied  di¬ 
rectly  to  these  artificial  time  series  to  obtain  estimates  of 
their  spectra  as  well  as  the  transfer  functions  for  each  chan¬ 
nel  and  desired  signals.  Since  all  the  ideal  spectra  are  known, 
the  effects  of  bias  errors  and  statistical  fluctuations  about 
the  mean  introduced  by  the  algorithm  can  be  explicitly  displayed 
and  studied  parametrically  not  only  in  terms  of  the  magnitude 
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of  M  but  also  as  a  function  of  specified  signal-to-noise  ratios 
Such  a  numerical  simulation  is  the  subject  of  this  paper. 

1.3  SCOPE  OF  PARAMETRIC  NUMERICAL  SIMULATIONS 

The  scope  of  these  parametric  investigations  is  limited  to 
preliminary  tests  of  the  algorithm  itself  under  arbitrary  but 
known  conditions,  and  it  deliberately  excludes  any  attempt  to 
evaluate  the  results  based  on  actual  data  gathered  by  the  Physi 
cal  Dynamics  experiments. 

The  subsequent  sections  of  this  paper  deal  with  the  follow 
ing:  (1)  a  brief  outline  of  the  mean-squared  estimation  algo¬ 

rithm  for  a  multiple-input  linear  system,  (2)  detailed  descrip¬ 
tion  of  the  numerical  simulation  scheme  for  a  three-channel 
linear  system,  (3)  presentation  of  the  numerical  results  and 
finally,  (4)  the  concluding  remarks. 
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2.  MEAN-SQUARED  ESTIMATION  ALGORITHM 


The  formulation  of  the  mean-squared  estimation  algorithm 
used  in  the  Physical  Dynamics  experiments,  as  described  in  the 
previous  section,  is  now  considered  in  order  to  bring  out  its 
basic  limitations.  This  section,  therefore,  considers: 

(a)  basic  concepts  in  the  formulation  of  the  algorithm,  (b)  re¬ 
construction  of  the  spectrum  for  the  main  data  channel,  (c)  fidel¬ 
ity  of  the  reconstructed  spectrum  of  the  imbedded  signal  in  the 
main  channel,  (d)  constraints  from  the  usual  assumptions  concern¬ 
ing  the  functional  linearity  among  the  superposed  time  series 
in  the  main  data  channel  and  those  from  the  subsidiary  measurements, 
(e)  constraints  from  the  frequency  averaging  of  auto-  and  cross- 
spectral  densities  as  used  in  common  practice,  and  (f)  formulation 
of  the  algorithm  for  a  subsidiary  nine-channel  system  as  used  in 
the  Physical  Dynamics  experiments. 

2.1  BASIC  CONCEPTS 

The  basic  concepts  underlying  the  reconstruction  algorithm 
are  illustrated  in  Pigs,  la,  b;  where  Fig.  la  (left)  shows  both 
the  main  data  channel  (top)  as  well  as  the  subsidiary  channels 
(bottom)  for  the  reconstruction  of  the  main  data  channel.  The 
time  series  y(t)  represents  the  main  data  channel,  which  is 
connected  to  the  primary  sensor  or  gradiometer  output  in  these 
particular  considerations.  The  main  data  y(t)  is  given  by 

q 

y (t)  *  E,  (t)  +  5Z  Mt)  +  z(t)  >  (12) 

1  i=2  1 
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MAIN  DATA  CHANNEL,  y(t) 


Schematic  representation  of  linear  system  utilized  in  the 
algorithm  to  reduce  the  Physical  Dynamics  measurements  in 
a  gradiometer  output,  (a)  Representation  of  superposed 
time  series  in  main  data  channel,  y(t),  and  time  series 
in  output  of  subsidiary  channels,  (b)  Reconstruction 
of  main  data  channel  for  simplified  functional  linear  re¬ 
lationship  (Eq.  15)  among  superposed  time  series  in  main 
data  channel  and  time  series  in  subsidiary  channels. 


i.e.,  by  the  desired  signal  denoted  by  the  time  series  ^(t) 
and  a  superposition  of  time  series  arising  from  several  but  not 
necessarily  independent  physical  sources  ^2 ) . . . S  (t ) ,  as  well 
as  a  composite  noise  process  z(t),  whose  constituent  time  series 
are  unknown.  Physically,  the  signal  £^(t)  represents  that  part 
of  the  voltage  fluctuations  generated  in  the  main  data  channel 
y(t)  due  to  the  interaction  of  ocean  currents  with  the  earth's 
magnetic  field.  Hence,  the  remaining  voltage  fluctuations  of 
y(t)  from  the  other  time  series  in  the  right-hand  side  of  Eq.  (12) 
constitute  sources  of  "noise"  with  respect  to  the  signal  £^(t). 

The  objective  of  the  reconstruction  algorithm  is  to  isolate  the 
signal  time  series  through  utilization  of  the  subsidiary  data 
channels  given  by 

x ( t )  =  [x1(t),  x2(t) . . .x1(t) . . . .x  (t)]  ,  (13) 

which  is  a  q-dimensional  input  vector  describing  the  subsidiary 
data  channels  indicated  in  Pig.  la.  Each  of  the  subsidiary  time 
series  in  Eq.  (13)  is  a  result  of  a  measurement  of  a  physical 
process  which  is  known  to  be  linearly  related  to  the  superposed 
time  series  £^(t),  )  •  •  •  •  (t ) . . .  .  £  (t )  in  the  main  data 

channel  y(t).  The  subsidiary  time  series  are  modeled  as  sta¬ 
tionary  stochastic  processes  with  statistical  properties  that 
need  not  go  beyond  the  second  moments  due  to  the  interest  in 
measurements  of  the  spectra  of  the  desired  signal.  If  the 
linear  relationship  among  the  subsidiary  and  superposed  time 
series  in  the  main  data  channel  is  time  invariant,  they  are, 
in  general  related  by 

q  rco 

£  (t)  =2  1  win<T)  xn(t-T)dT,  i=l,  2...q  (14) 

1  n=l 

where  (t)  represents  the  impulse  response  function  in  the 
presence  of  the  noise  z(t)  and  Eq.  (14)  is  the  linear  super¬ 
position  of  the  convolution  integrals.  Physically,  Eq.  ( 1 4 ) 
indicates  that  each  of  the  superposed  time  series  (£^)  in  the 
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main  data  channel  is  the  results  of  a  linear  superposition  of 

the  contributions  to  it  from  each  of  the  q  subsidiary  data 

q 

channels  (  given  by  future  and  past  inputs  (-»  <  T  <  «) 

from  each  subsidiary  time  series.  Hence,  the  outputs  of  the 

subsidiary  and  main  channels  can  be  used  to  estimate  £;  (t) . 

£  (t),  provided  the  q  x  q  matrix  of  impulse  reponse  functions 
v*!  (t)  can  be  determined  from  the  experimental  data.  However, 
following  common  practice  (e.g..  Ref.  2),  the  algorithm  used  in 
the  Physical  Dynamics  interpretation  of  the  data  utilizes  a 
special  form  of  Eq.  (14),  viz.. 


€±(t ) 


V±(l)  X^t-TjdT, 


1*1,  2 . . . q  (15) 


wherein  only  q  impulse  response  functions  w^(t)  must  be  deter¬ 
mined  from  the  experimental  data.  Unlike  the  more  general 
linear  relationship,  Eq.  (14),  the  special  form  given  by  Eq.  (15) 
implies  a  one-to-one  correspondence  between  each  x^ ( t )  subsidi¬ 
ary  channel  and  its  corresponding  ?1(t)  time  series  in  the  main 
channel.  Physically,  the  specialized  linear  assumption  given 
by  Eq.  (15)  implies  two  conditions  (Ref.  2):  (a)  that  the 

output  of  several  simultaneous  inputs  to  the  main  channel  is 
equal  to  the  sum  of  outputs  produced  by  each  input  alone  and 
(b)  that  the  output  of  the  product  of  a  constant  times  an  in¬ 
put  is  equal  to  the  constant  times  the  output  of  the  input 
alone.  Hence,  Eq.  (15)  assumes  a  complete  functional  decoupling 
between  £^(t)  and  x  (t)  whenever  n  ^  i,  as  shown  in  Fig.  lb. 
Equivalently,  the  impulse  response  functions  in  Eq.  (14)  become 
related  to  those  in  Eq.  (15)  by 


w 


w1(t)  ;  i  *  n 
0  ;  i  ¥■  n 


(16) 


i.e.,  when  only  the  terms  along  the  principal  diagonal  of  the 
q  x  q  matrix  of  w^n(x)  are  non-zero.  Note  that  this  functional 
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decoupling  does  not  imply  statistical  independence  among  the 
various  time  series;  i.e.,  xn  and  are  not  necessarily  sta¬ 
tistically  independent  even  though  Eq.  (15)  holds.  Whether 
the  subsidiary  channels  are  related  to  the  constituent  compon¬ 
ents  in  the  main  channel  by  Eq.  (15)  or  by  the  more  general 
form  of  Eq.  (14)  depends  on  a  priori  knowledge  of  the  physical 
processes  modeled  by  each  x  (t),  and  is  not  readily  resolvable 
on  the  basis  of  the  experimental  data  alone.  Some  consequences 
of  assuming  incorrectly  that  Eq .  (15)  holds  always  will  be 
discussed  later.  The  fidelity  of  reconstruction  of  the  main 
channel  output  y(t)  in  terms  of  the  subsidiary  data  records 
xfi(t)  does  not  depend  on  the  choice  between  Eqs.  (14)  and  (15) , 
a  choice  that  becomes  important  only  in  the  subsequent  identi¬ 
fication  of  £^(t)  as  t*ie  desired  reconstructed  signal.  While 
Eq.  (15)  would  isolate  the  signal  £^(t),  a  reconstruction  based 
on  Eq.  (14)  would  yield  an  unknown  linear  superposition  of  the 
signal  and  "noise"  outputs  C?(t) . £  (t)  in  Eq> 

Finally,  it  must  be  emphasized  that  an  alternative  theoret¬ 
ical  procedure  to  obtain  the  desired  signal  (Fig.  1)  would 
be  to  (a)  utilize  only  the  subsidiary  data  channels  [e.g., 
x^(t)]  simulating  internal  wave  phenomena  and  (b)  use  the  Max¬ 
well  equations  to  calculate  directly  the  transfer  functions 
H1(f),  i.e.,  to  calculate  the  magnetic  gradients  induced  by 
internal  waves.  However,  this  more  complicated  approach  would 
utilize  only  the  available  multi-channel  hydrodynamic  data  in 
the  Physical  Dynamics  experiments  and  would  not  represent  an 
experimental  verification  of  the  gradiometer  measurements. 

The  latter  consideration  is  important  because  the  calculation 
of  this  transfer  function  must  also  utilize  auxiliary  assump¬ 
tions  concerning  the  distribution  of  internal  wave  energy  among 
wave  modes  as  well  as  the  directional  properties  (isotropic 
versus  nonisotropic)  of  the  relevant  processes  in  a  shallow 
ocean. 
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2.2  RECONSTRUCTION  OF  SPECTRUM  FOR  MAIN  DATA  CHANNEL 


Based  on  Eq.  (15),  Fig.  lb  illustrates  the  reconstruction 

A 

or  estimation,  y(t),  of  the  time  series  y(t)  in  the  main  data 
channel  for  z(t)  =  0.  If  w^t)  =  hi(r)  for  this  latter  condi¬ 
tion,  the  estimator  y(t)  is  given  by  Eqs.  (12)  and  (15)  as 

q  roc 

y(t)  =  ^  I  h. (t)  x. (t-T)di  (17) 

i  =  l  J-co  1  1 

i.e.,  each  subsidiary  channel  x^(t)  is  passed  through  a  time 
invariant  linear  filter  with  impulse  response  function  h. (r) 
and  the  results  summed  to  yield  the  estimate  y(t).  The  criterion 
for  good  reconstruction  of  y(t)  is  taken  to  be  the  statistical 
average  of  the  mean-squared  error  or  difference  between  the 

A 

measured  time  series  y(t)  and  the  estimator  y(t),  i.e., 

<e2>  =<  |  y  ( t )  -  y (t )  | 2  >  .  (18) 


The  correct  impulse  response  functions  hj_(t)  in  Eq.  (17)  are 
chosen  such  that  <  e2  >  is  minimized.  If  the  stochastic  processes 
are  stationary,  it  is  found  that  the  hi(x)  functions  that  yield 
such  a  minimum  must  satisfy  the  system  of  integral  equations 


hj  (t'  )R^1(t-t'  )dT ' 


(t). 


1  =  1,2 


q  (19) 


where  and  Ryi  are  auto-  and  cross-correlation  functions  de¬ 
fined  by  Eqs.  (2)  and  (4),  i.e.. 


R11(t)  *  <  x1(t+x)  x±(t)  > 
Rji(f)  B  <  Xj ( t+T )  x±(t)  > 
Ryl(x)  =  <  y(t+-r)  xA(t)  > 


The  q  impulse  response  functions  satisfying  Eqs.  (19)  de¬ 
fine  a  q-dimensional  Wiener  filter.  Since  all  physically  real¬ 
isable  impulse  response  functions  must  vanish  for  the  negative 
values  of  their  argument  x,  the  limits  of  integration  must  be 
truncated  to  run  from  0  to  As  a  result,  the  matrix  integral 
equation  is  then  of  the  Wiener-Hopf  type;-  the  solution  of  which, 
in  general,  presents  difficult  problems.  However,  if  one  is 
willing  to  accept  impulse  response  functions  which  are  non- 
causal,*  i.e.,  h^  (t)  /  0,  x  <  0,  the  solution  of  Eq.  (19)  is 
found  by  a  straightforward  application  of  Fourier  transforms. 
Thus,  in  the  frequency  domain,  Eq.  (19)  factors  into 


q 

53  H.(f)  S1t(f) 

J-l  J  J 


syi(f) 


i  =  1,  2 _ q  (20) 


where  Hj(f),  S^(f)  and  S  (f)  are  Fourier  transforms  of  h^  (t), 
R  (t)  and  Ryi<x),  respectively.  Fortunately,  the  possible 
noncausal  nature  of  hj  (x)  as  obtained  from  Eq .  (20 )  does  not 
appear  to  be  too  important  in  the  digital  realization  of  the 
filter  functions,  so  that  there  is  no  need  to  deal  with  the 
substantially  more  complicated  Wiener-Hopf  solution. 

In  order  to  expand  Eq.  (20)  in  matrix  form,  it  becomes 
convenient  to  introduce  the  following  notation:** 


H(f)  =  [Hx<f ) ,  H2(f)  ... 

Vf>  -  tvf><  Vf> 


Vf)] 

(21) 

Qy(f)] 

(22) 

A  filter  is  said  to  be  noncausal  if  it  can  produce  an  output 
before  an  input  is  applied. 

ft  ft 

Since  the  algorithm  is  described  in  some  detail  in  Ref.  2  ,  the 
same  notation  is  used  here  as  far  as  possible  to  facilitate 
the  interpretation  of  comparable  results. 


15 


and 


Sxx<f>  - 


n(f> 

s12(f)  .... 

••••  Slq ( f ) 

21(1) 

S22 ( f )  •••• 

....  s2  (f) 

qi(f) 

s  2(f)  .... 

....  sqq(f) 

(23) 


where  H(f)  Is  the  q-dimensional  frequency  response  or  transfer 

function  vector  corresponding  to  that  in  Eq .  (13)  under  the 

assumption  of  Eq.  (15),  3  the  q-dimensional  cross-spectrum 

xy 

vector  of  the  subsidiary  inputs  x^(t)  with  the  main  channel  time 
series  y(t),  and  §  (f)  the  q  x  q  spectral  matrix  of  the  subsi¬ 

diary  time  series.  The  notation  in  the  right-hand  side  of  Eq.  (22) 

+  -+■ 

has  been  simplified  by  replacing  Sx  y  with  where  i  =  1,  2  ...q. 

Similarly,  in  Eq .  (23)  SY  v  is  replaced  with  S. . ,  where 

i  J  ^<3 

i,  j=l,  2....q.  It  should  be  noted  from  Eqs.  (6)  and  (8)*  that 
equation  (23)  is  a  Hermitian  matrix  [S^j ]  characterized  by  real¬ 
valued  elements  along  the  principal  diagonal  (i=J)  and  by  complex 
conjugates  for  any  two  elements  (S^j  ,  Sj^)  situated  symmetrically 
with  respect  to  the  principal  diagonal.  A  prime  superscript  will  be 
used  to  denote  the  transpose  of  a  given  vector;  i.e.  H'  to  indi¬ 
cate,  for  example,  the  vertical  column  corresponding  to  the  row 
given  by  Eq .  (21).  Finally,  a  t  superscript  will  indicate  the 
complex  conjugate  transpose  matrix  (adjoint)  to  denote  from  Eq. 

"►  i 

(22),  for  example,  by  Si  the  vertical  column  with  elements 

xy 

S*iy...S*qyj  whereas  an  asterisk  superscript  will  denote,  as  be¬ 


fore,  the  conjugate  of  a  complex  variable.** 
preceding  notation,  Eq.  (20)  reads 


In  terms  of  the 


H  S 


xx 


=  S* 


xy 


(24) 


~K 

These  equations  can,  of  course,  be  generalized  to  include  any 
Xj_  in  Eq.  (6)  and  any  combination  of  X*j_Xj  for  i  /  J  in  Eq.  (8). 

To  simplify  the  already  cumbersome  notation,  the  explicit  de¬ 
pendence  on  frequency  (f)  will  be  omitted. 
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i 


from  which  the  solution  for  the  transfer  function  column  matrix 
is 


H1 


=  (S  ) _1  s' 
x  xx '  xy 


(25) 


where  the  symmetries  S., 


=  S*  have  been  used. 


,  =  S#  and  S 
ij  ji  yn  ny 

The  filter  transfer  functions  given  by  Eq.  (25)  yield  the 


smallest  mean-squared  errors  in  Eq.  (19) 
error  is  given  by 

Q 


This  mean-squared 


<e2> 


min 


/•OO 

■/-  -  s 


sJy  Vdf 


(26) 


The  integrand  in  this  expression  is  just  the  spectrum  of  the 
residual  process  (uncancelled  noise)  y(t)  -  y(t)  in  the  main 
data  channel.  Upon  substituting  the  solution  from  Eq.  (25) 
for  H*  into  Eq.  (26),  one  can  rewrite  the  latter  entirely  in 
terms  of  the  auto-  and  cross-spectral  density  functions,  i.e.. 


<e 


/■•  r  S*  (S  )_,s'  1 

2>  =/  s  i  -  '-** - audf 

min  yy  L  syy  J 


or 


<e2> 


min 


-r 


s  [1  -  Y2]  df 
yy 


(27) 


where  the  quantity  y2  is  known  as  the  multiple  coherence  func¬ 
tion  as  given  by 


S*  (S  )"lS' 

'  w  '  yy 


y2  =  xy  xx 
Y  S 

yy 


(28) 


and  is  equal  to  the  normalized  cross  spectral  density  between 
the  main  channel  time  series  y(t)  and  the  estimator  y(t): 
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(29) 


y2  = 


is*  I2 
_LizL. 
s  s** 
yy  yy 


If  the  main  data  channel  time  series  y(t)  were  completely 
decomposable  into  a  sum  such  as  Eq.  (17),  then  one  would  neces¬ 
sarily  have  y2  =  1  and  the  mean  squared  error  would  be  zero. 
Conversely,  the  condition  <e2>mln  -  0  yields  from  Eq.  (27) 

Y2  =  1,  because  (excluding  the  trivial  case  of  S  =0) 

y  y 


and  1  - 
that 


Y2  are  always  positive, 


yy 


yy 


=  s 


xy 


(s  )-1s 

xx  xy 


From  Eq.  (28),  it  would  follow 

(30) 


which,  with  Eqs.  (24)  and  (25)  yield 


yy 


s  H  Sxx  H 


(31) 


or 


S11 

S12  • 

. ...  slq 

S21 

S22  * 

. ...  s2q 

H* 

(31a) 

syy 

-  [«!••• V 

• 

• 

• 

SQ1 

8Q2  ' 

*  •  •  •  s 

qq 

H* 

q 

By  direct 

calculation 

,  it 

can  be 

shown  that  the  right-hand  side 

of  Eq.  (31)  is  just  the  spectrum  of  the  estimator  y(t),  i.e., 

s0v"  Hence  y2  =  1  implies  that  i.e.,  the  spectra  of 

jjf  yy 

the  time  series  y(t)  and  the  estimator  y(t)  are  identical. 

Thus,  the  departure  of  y2  from  unity  may  be  taken  as  a  measure 

of  the  degree  to  which  the  spectrum  of  the  estimator  fails  to 

approximate  the  spectrum  of  the  time  series  y(t)  in  the  main 

channel . 
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2.3  FIDELITY  OF  THE  RECONSTRUCTED  SPECTRUM  OF  THE  IMBEDDED  SIGNAL 

Returning  to  Eq.  (12),  the  main  problem  of  interest  is 
the  decomposition  of  the  series  y(t)  into  its  constituents  so 
as  to  estimate  the  spectrum  of  the  signal  ?^(t).  Clearly, 
unless  one  can  associate  particular  subsidiary  channels  with 
the  correct  set  of  superposed  series  £  (t)  in  the  main  data 
channel,  no  such  decomposition  is  possible.  Suppose  the  assump¬ 
tion  is  made  that  each  E^(t)  is  uniquely  related  to  ,  as  in 
Eq.  (15).  Obviously,  if  in  addition  z(t)  =  0  in  Eq.  (12),  then 
w^(t)  =  hi(t),  which  are  just  the  estimated  impulse  response 
functions  and  whose  Fourier  transforms  are  given  by  Eq.  (25). 
Hence, 

»  f  h± ( t— t )  xi(x)dT  .  (32) 

—00 

In  particular,  the  reconstructed  signal  spectrum  is 


S 


«1«1 


Hi  I2 


’ll 


(33) 


When  y(t)  comprises  a  noise  output  z(t)  not  accounted  for  by 
the  x1(t)  subsidiary  series,  the  use  of  Eq.  (33)  will  involve 
errors  since  perfect  reconstruction  of  y(t)  would  no  longer  be 
possible.  While  the  closeness  of  the  multiple  coherence  func¬ 
tion  to  unity  may  be  regarded  as  a  figure  of  merit  of  the 
fidelity  of  reconstruction  of  the  overall  output  y(t),  it  is 
not  necessarily  a  good  measure  of  the  fidelity  of  reconstruc¬ 
tion  of  a  particular  constituent  of  y(t),  e.g.,  E^(t)  in  Eq. 
(12).  The  quantity  that  plays  a  role  analogous  to  the  multiple 
coherence  for  Syy  and  gives  a  quantitative  measure  of  the 
fidelity  of  reconstruction  of  is  the  so-called  partial  (or 
conditional)  coherence  function.  This  function  is  defined  as 
the  normalized  cross-spectral  density  between  a  pair  of  condi¬ 
tioned  random  variables  via  the  following  procedure. 
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Instead  of  estimating  y(t)  on  the  basis  of  the  q-dimen- 
sional  vector  given  by  Eq.  (13),  the  procedure  consists  in 
doing  so  only  on  the  basis  of  q  -  1  subsidiary  channels. 
Specifically,  if  the  excluded  channel  output  is  x.(t),  the 
estimator  y  (t)  based  on  the  remaining  channels  is  given  by 

y(l)(t)  =  i  f  ( x )  x.  (t  -  x)  dx  .  (34) 

1=2  yl  1 

The  impulse  response  functions  h^j  (t)  are  now  determined 
subject  to  the  minimization  of 

<e(l)%  =  <  | y ( t )  -  y(l)(t)|2  >  .  (35) 

If  the  Fourier  transform  of  hy.^  (x)  is  denoted  by  Hy^  and  the 
corresponding  column  matrix  by  the  solution  to  the  mini- 

mization  problem  is  again  given  by  a  matrix  equation  similar  to 
Eq.  (25),  viz. 

•*( ,  ^ 

where  S' '  is  the  reduced  (q  -  1)  x  (q  -  1)  covariance  matrix 
obtained  from  S  by  deleting  the  first  row  and  first  column. 

■*■.(  i 

Similarly,  S A.  ’  is  the  reduced  column  matrix  with  q  -  1 
elements  corresponding  to  the  horizontal  vector: 


sqy^ 


(37) 


Having  determined  y^(t),  the  reduced  (conditioned)  random 
process  A^y(t)  may  be  considered  by  subtracting  from  y(t)  its 
best  linear  estimate  based  on  the  knowledge  of  XgCt),  x^(t) 

.  Xq (t),  1 . e . , 


Axy(t)  *  y(t)  -  y^(t)  . 


(38) 
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Clearly,  with  the  best  linear  estimate  obtained  from  X£(t) 
.  xq(t)  subtracted  out,  ^yCt)  should  be  a  better  repre¬ 
sentation  of  §^(t)  than  y(t)  itself  provided,  of  course,  that 
y(t)  indeed  does  contain  a  linear  functional  dependence  on 
x^(t),  i.e.,  that  the  assumption  given  by  Eq.  (15)  holds.  In 
the  same  manner  as  A^y(t),  one  can  also  estimate  x^(t)  on  the 

basis  of  x^(t)  .  x^(t)  and  subtract  this  estimate  from 

x^t),  thereby  forming  the  reduced  process  A^x(t),  i.e.. 


A-j^xCt)  =  x-^t)  -  x1(t) 


where 


4  oo 

c-^t)  =  S  f  (t)  X^t-OdT 

1=2  J-co  x 


If  Fourier  transforms  of  the  impulse  response  function  hr  '  (t) 
/  . \  x,  1 


are  denoted  by  H'  i  with  the  corresponding 

-*■,  ( i )  xi  -1 


column  matrix 


,  one  has 


5;(°  -  feV  w 


where  Sx^  '  is  the  reduced  column  matrix  with  q  -  1  elements 
corresponding  to  the  horizontal  vector: 


Sxl’  ■  rs21  s31 


Since  Ax^t)  is  formed  by  subtracting  from  x1(t)  the  best 

linear  estimate  x^(t)  from  Xg(t)  .  x^Ct),  it  is  reasonable 

to  (a)  associate  a  high  degree  of  correlation  between  the 
reduced  (conditioned)  time  series  A1y(t)  and  A1x(t),  and  (b) 
expect  a  faithful  reproduction  of  the  signal  spectrum  as  given, 
for  example,  by  Eq.  (33).  The  correlation  between  the  residual 
A1y(t)  and  A^(t)  series  can  be  expressed  in  the  frequency 
domain  as  a  partial  coherence  function.  The  auto-  and 


cross-spectral  densities  of  the  residual  processes  A-^xtt) 
A.y(t)  are  defined  by 

•1*  CD 


syy;  2,3 . q  *  fj'W 

00 

Sll;  2,  3,  .  q  ~  ^^[<A1x(t+T)A1x(t)>]e 

00 

Sly ;  2,  3,  .  q  =  f_J~<L ix(t+T)Aiy(t)>^e 

The  partial  coherence  function  is  then  defined  by 

Y2  _ \siy;2,3, . g|2 _ 

ly;  2’3,  .  q  syy;  2,3 . .  sii;  2,3,. 


dx 


dx 


dx 


and 


( ^3) 

(44) 

(45) 


_  (46) 

q 


Successive  substitution  for  Eqs.  (38),  (39),  (34),  (40),  and  (36) 
for  A1y(t),  A^xft),  y^(t),  x-^t)  and  Hy+^  \  respectively,  into 
the  definitions  (43)  -  (45)  yield  for  the  latter  the  following: 


syy;  2,3, 
snj  2,3, 
siy;  2,3, 


In  the  forms  given  by  Eqs.  (47)  -  (49),  the  conditional  spectral 
densities  and  hence  the  partial  coherences  may  be  computed 
directly  from  the  unconditional  spectral  densities  of  the  x(t) 
and  y(t)  time  series. 


The  conditional  spectral  densities  and  partial  coherences 
may,  of  course,  be  defined  with  respect  to  any  one  of  the  sub¬ 
sidiary  q-channels  in  Eq.  (13).  A  convenient  way  of  keeping 
track  of  the  indexes  in  the  corresponding  Eqs.  (47)  -  (49)  is 
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to  associate  the  first  row  and  column  of  the  subsidiary  matrix 
S  with  the  channel  of  interest  and  perform  a  cyclical  permu- 

A  A  2 

tation  of  the  remaining  indexes.  For  example,  Y0...  ,  >.  ■, 

refers  to  the  partial  coherence  relative  to  x2(t).  The  corres¬ 
ponding  conditional  spectral  densities  are  then  S 
etc.  Since  in  this  manner  all  the  formulas  involving  conditional 
spectral  densities  pertaining  to  the  channel  x^t)  can  readily 

be  modified  for  any  one  of  the  x2(t)  .  xq(t)  channels,  one 

can  continue  for  convenience  to  refer  specifically  only  to  the 
channel  x^t). 

Although  the  channel  transfer  functions  may  be  computed 
directly  from  Eq.  (25),  it  is  also  possible  to  compute  them 
from  the  conditional  spectral  densities.  The  required  formulas 
can  be  obtained  by  combining  Eq.  (25)  with  EqB.  (48)  and  (49). 

For  this  purpose,  the  matrix  relation  given  by  Eq.  (24)  is  first 
written  in  the  following  partitioned  form 


■  1 
sn  i 

S(o 

Slx 

!<■> 

XX 

H, 


V-1 


ly 

£»(*) 

xy 


(50) 


where  Is  the  q-1  dimensional  column  matrix  containing  the 

transfer  functions  [H^,  .  H*].  From  the  matrix  Eq.  (50! 

it  is  obtained 

4/n  . 

(51a) 
(51b) 

Eliminating  yields 


su  Hi* + 


o’  (O  H  »  +  S^l)  H  + 

Slx  H1  +  Sxx  Hq-1 


iy 


5»(0 

xy 


sly  -  53°  (^))“  «/. 


In  view  of  Eqs.  (48)  and  (49)  one  has 


*  _  ly;  2, 


s 


ll;  2,  3. 


or  generalizing  to  an  arbitrary  channel  X^t) 
IT  *  _  Siy;  i+l.i+2,  . q.1,2... 


*i  S 


ii;  i+l.i+2. 


Q  >  1  >  2, . i-1 


For  a  "perfect"  reconstruction  of  y(t),  the  partial  coher¬ 
ences  must  be  unity,  a  result  that  is  shown  by  assuming  first 
that  the  mean  squared  error  in  Eq.  (26)  is  zero.  Since  the 
integrand  is  positive,  it  follows  that 

Q 

Syy  *£  Sly  H1  (55> 


■i +s 


siy  Hi  ■  Hlsly  +  \-l  Sxy° 


On  the  other  hand,  from  Eq.  (51b)  it  is  found  that 

Vl  "  Sxy  \Sxx/  -  Sxl  Vsxx  )  H1  •  (57) 

The  combination  of  Eqs .  (56)  and  (57)  yields 

■yy-^Mr  (^,)',^')•(58, 


Substitution  from  Eqs.  (47)  and  (48)  gives  the  relationship 


syy;  2,  3, 

. q 

Hi  siy;  2, 

3, 

•  •  •  *  ' 

. .  q  . 

(59) 

view  of  Eq.  (53) 

there  also 

follows 

syy;  2,  3, 

.  q 

IhiT  su. 

2, 

3,  . 

(60) 

|siy;  2,  3,  .. 

....  qi*  • 

M!  4. 

2, 

3,  • 

•  •  •  ■  •  Q 

(61) 

so  that  the  partial  coherence  as  given  by  Eq.  (46)  is  identi¬ 
cally  unity.  Since  this  is  also  true  for  each  of  the  remaining 
q-1  subsidiary  channels,  it  is  seen  that  perfect  reconstruction 
of  y(t)  as  given  by  Eq.  (12)  implies  that  for  z(t)  0  all  the 
partial  coherences  must  be  unity.  It  is  therefore  common 
practice  to  associate  a  high  value  of  the  partial  coherence  for 
a  particular  channel  with  high  fidelity  of  spectral  reconstruc¬ 
tion  of  the  corresponding  signal  as  given,  for  example,  by 
Eq.  (33). 

2.4  CONSTRAINTS  FROM  USUAL  LINEARITY  FORM 

It  should  be  noted  that  the  foregoing  association  between 
fidelity  of  reconstruction  of  A^*)y(t)  and  high  values  of  the 
partial  coherences  Is  justified  only  If  the  linear  relationship 
(15)  is  known  to  hold  between  the  superposed  and  subsidiary  time 
series  £^(t)  and  x^(t),  respectively.  To  illustrate  this  point, 
consider  the  ideal  case  when  z(t)  =  0  in  Eq.  (12)  and  assume 
that  the  general  linear  relationship  given  by  Eq.  (14)  holds 
instead  of  Eq.  (15).  The  reconstruction  of  y(t)  is  then  given 
by 

q  q  00 

y(t)  =  £  J  win  <*)  xn  (t-T)  dT 
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or  putting  this  result  in  the  form  of  (17) 

q 


y<t)  ■  t-1  L  [n?!  t*(T>]  xl(t- 


t)  dx  . 


(62) 


The  algorithm,  based  on  Eq.  (15) »  will  then  estimate  the 
Fourier  transform  of  the  impulse  response  functions 

A 

(63) 


h.  (x)  =  iC  w  (x)  . 
i  n=l  ni 


Denoting  the  Fourier  transform  of  wni(t)  by  Hni(f),  then  the 
estimated  transfer  function  would  be 

q 

<1  ’  Hni  •  (61|> 


H. 


n«l 


Thus,  if  the  signal  were  associated  with  the  first  subsidiary 
channel,  Eq.  (33)  would  yield  the  following  estimated  spectrum 
of  the  signal 

q 

£  « 


£i«i 


n=l 


nl 


'11 


(65) 


However,  the  actual  signal  imbedded  in  y(t)  is,  from  (14), 

q 


(t )  ]actUai 


=  ^  /  w  (x)  x.  (t-x)dx 

i=l  11  1  ' 


(66) 


with  the  corresponding  spectrum 

q  q 


^1^1^  '*actual 


=  £  £ 

i=l  n=l 


Hli  Hln  Sin 


(67) 


Clearly,  Eq.  (67)  can  be  substantially  different  from  Eq.  (65). 
It  should  be  noted  that  under  the  stipulated  conditions,  the 
multiple  coherences  and  each  of  the  partial  coherences  are 
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necessarily  equal  to  unity.  Hence,  in  absence  of  a  priori 
information  on  the  interchannel  functional  linear  relationships, 
the  conclusion  would  be  made  with  "great  confidence"  that 
Eq.  (65)  is  indeed  the  signal.  With  the  algorithm  in  its 
present  form,  no  general  remedy  is  apparent  to  remove  this 
difficulty.  In  essence,  the  degree  of  confidence  assigned  to 
the  correctness  of  the  reconstructed  signal  can  then  be  no 
greater  than  the  degree  of  confidence  in  the  a  priori  knowledge 
that  a  subsidiary  channel  i  responds  only  to  one  physical 
mechanism.  These  considerations  lead  then  to  a  need  for  a 
linear  superposition  into  one  channel  of  all  the  subsidiary 
time  series  that  are  functionally  related  to  the  signal,  and 
that  may  be  suspected  to  be  present  in  more  than  one  subsidiary 
channel.  This  linear  superposition  would  then  yield  a  smaller 
set  of  subsidiary  channels  that  are  indeed  functionally 
decoupled.  One  way  to  test  the  functional  relationship  among 
the  subsidiary  channels  would  obviously  be  to  test  the  assump¬ 
tion  implied  by  Eq.  (15)  with  the  outputs  of  the  subsidiary 
channels  by  themselves. 

2.5  FORMULATION  OF  NINE-CHANNEL  LINEAR  SYSTEM 

The  formulation  of  the  reconstruction  algorithm  for  the 
nine  subsidiary  channels  of  the  Physical  Dynamics  experiments 
is  based  on  Eq.  (15),  Instead  of  Eq .  (14),  for  the  functional 
linear  relation  among  the  superposed  time  series  in  the  main  data 
channel  [e.g.,  Eq.  (12)]  and  the  subsidiary  channels  [Eq.  (13)3  - 
In  view  of  the  typical  result  given  by  Eq.  (33) j  which  relates 
the  spectrum  of  a  given  superposed  time  series  in  the  main 
channel  to  the  transfer  function  of  the  corresponding  subsidiary 
channel,  the  main  interest  in  this  formulation  is  focused  on 
the  transfer  functions  of  the  subsidiary  channels.  Because  the 
fidelity  of  reconstruction  of  the  superposed  time  series  in  the 
main  data  channel  is  given  by  the  corresponding  partial  coher¬ 
ences  [e.g.,  Eq.  (46)3,  the  considerations  on  the  transfer 
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functions  for  the  subsidiary  channels  must  include  the  corre¬ 
sponding  partial  coherences.  Hence,  the  formulation  of  the  nine- 
channel  linear  system  is  based  on  the  specialization  of  Eq.  (54) 
for  the  transfer  function  for  q  =  9  as  well  as  the  corresponding 
generalized  Eq.  (46)  for  the  corresponding  partial  coherence. 


The  formulation  of  the  nine-channel  linear  system,  based 
on  the  results  of  Ref.  3,  utilizes  Eqs .  (13),  (21),  (22)  and 
(23)  for  the  vectors  x(t),  i?(f),  §  and  2  ,  respectively. 

.  xy  XX 

The  q-dimensional  vector  S  consists  of  cross-spectra  densities 

xy 

involving  each  of  the  subsidiary  channels  with  the  main  data 
channel,  whereas  the  §  is  a  q  x  q  Hermitian  matrix  involving 

X  X 

only  the  subsidiary  channels.  As  indicated  in  Refs.  2  and  3, 

the  use  of  Eq.  (54)  for  the  transfer  function  as  well  as  of 

Eq.  (46)  for  the  partial  coherences  may  be  put  in  a  partitioned 

form  similar  to  that  of  the  matrix  Eq.  (50).  This  procedure 

requires  the  partitioning  of  an  augmented  matrix  2  (f) 

yxx 

given  by 
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The  partitioned  matrix  (68 ) 

can 

be 

represented  as 
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where  matrices  A,  B,  C  and  D  with  their  specified  ranks  corre¬ 


spond  to  each  respective  matrix  in  the  partition  of  (68).  It 

is  to  be  noted  that  matrix  D  corresponds  to  matrices  C 

■*'  (1)  xx 

and  B  are  augmented  matrices  involving  &Xy  »  and  matrix  A  is 

augmented  by  S  ,  Syj_  and  .  The  partitioning  of  Eq.  (69) 
permits  the  definition  of  a  residual  matrix  2  #  -  _  Q  given 

by  (Refs.  3,  4): 


Ky,2,3;..9  *  A  -  BD“C  <7°> 

where  D_1  is  the  inverse  matrix  of  D.  The  residual  matrix 
3  . 2  -3  a  is  the  same  rank  as  matrix  A  (2x2),  with 

elements  given  by 
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xy ;  2 , 3 . 


•  9 


syy;2,3. .  .9 


syi;2,3. 


iy;2,3. ♦  -9 


11;2,3. . .9 


(71 ) 


i.e.,  by  the  elements  given  by  Eqs .  (47)  to  (49).  Therefore, 
it  should  be  noted  that  the  second  row  of  (71 )  contains  the 
numerator  and  denominator  for  the  transfer  function  as  given  by 
Eq.  (53)  >  whereas  the  two  diagonals  contain  the  denominator 
(principal  diagonal)  and  numerator  of  the.  partial  coherence  as 
given  by  Eq.  (46).  Hence, 


and 
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(72) 


Y  ly;2,3. . .9  *  TS 


■.-Siyj.^.---9I_2 _ 

ii;2,3. . .9)(syy;2,3. . .9 


T 


(73) 


29 
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The  subsequent  transfer  functions  and  partial  coherences 

.  Y,9y;l,2...8  are  emulated  successive 

permutations  of  the  subindexes  in  the  augmented  matrix  (68 ) . 

The  augmented  matrix  (68)  contains  all  the  spectral  density 
functions  required  for  the  calculation  of  the  transfer  functions 
and,  hence,  the  corresponding  spectra  of  the  superposed  series 
in  the  main  data  channel  (e.g.,  Eq.  33).  These  spectral  density 
functions  are  in  turn  determined  from  the  corresponding  esti¬ 
mators  as  given  by  Eqs .  (5)  to  (11),  a  procedure  that  requires 
use  of  the  output  time  series  from  the  respective  subsidiary 
[Eq.  (13)]  and  main  data  channels. 


3.  SIMULATION  OF  A  THREE- CHANNEL  LINEAR  SYSTEM 


f 


A  test  of  the  algorithm  utilized  for  the  reconstruction  of 
the  main  data  channel  in  the  Physical  Dynamics  experiments  (Pig.  1) 
can  be  performed  by  simulating  the  data  records  for  the  main  data 
channel  and  only  two,  instead  of  the  nine,  subsidiary  channels 
used  in  such  experiments.  By  expanding  such  approach  to  include 
a  successively  increasing  number  of  subsidiary  channels,  it 
would  become  possible  to  identify  any  systematic  growth  in  the 
bias  of  the  spectrum  of  the  imbedded  signal  as  a  result  of  (a) 
the  increased  number  of  estimators  in  the  augmented  matrix  in 
Eq.  (68)  and  (b)  the  algebraic  operations  [Eqs.  (72  ),  (73  )]  with 
such  biased  estimators  during  the  application  of  the  algorithm 
to  a  larger  number  of  channels. 

The  test  of  the  algorithm  based  on  the  use  of  the  main  data 
channel  [y(t)]  and  two  subsidiary  channels  [x^t)),  x2(t)]  involves 
the  construction  of  two  3x3  augmented  matrices  [Eq.  (68)]  as 
follows:  (a)  one  characterized  by  ideal  spectra  and  zero  cross¬ 

spectra  densities  (S12  =  S21  =  0)  for  the  uncorrelated  subsidiary 
series  x-^(t)  and  x2(t).  This  matrix  yields  then  the  ideal  trans¬ 
fer  function  [Eq.  (72)]  and  other  parameters  such  as  the  ideal 
coherence  based  on  the  use  of  ideal  spectra  [e.g.  Eq.  (73)], 
and  (2)  another  one  characterized  by  estimators  [Eq.  (11)]  with 
non-zero  cross-spectra  for  the  subsidiary  channels,  depending 
upon  the  degree  of  frequency  averaging.  This  augmented  matrix, 
based  on  estimators,  yields  the  estimated  transfer  functions 

A  A  A 

(Pig.  1)  and  the  reconstructed  signal  spectrum  =  |  | 2 

[Eq.  (33)].  The  ultimate  test  of  the  algorithm  is  then  given  by 
comparing  the  fidelity  of  the  reconstructed  signal  spectrum  as 
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well  as  the  partial  coherence,  i.e.,  the  criterium  for  the  fidel- 

A 

ity  of  the  reconstructed  time  series  [y(t)]  in  the  main  data 
channel  (Fig.  1).  A  more  detailed  description  of  the  construc¬ 
tion  of  these  augmented  matrices  as  well  as  the  specification 
of  the  ideal  signal  spectrum  in  the  gradiometer  output  is  given 
in  the  sequel. 


The  construction  of  the  augmented  matrices  require  deter¬ 
mination  of  the  corresponding  spectral  density  functions  from 
time  series  for  the  gradiometer  output  and  subsidiary  channels 
[Eqs.  (5)  -  (11)].  The  simulation  of  data  records  for  the  gra¬ 
diometer  [y(t)]  and  two  subsidiary  channels  [x-^t),  x2(t)]  is 
illustrated  schematically  in  Fig.  2,  which  also  shows  four 
filters  denoted  as  G^f),  G2(f),  K^Cf)  and  K2(f).  Each  filter 
is  characterized  by  its  bandwidth  and  the  center  frequency  for 
peak  amplitude  of  the  spectral  density.  The  first  two  filters 
are  for  a  two-input  gradiometer  and  the  latter  two  for  the  subsi¬ 
diary  single-input  channels  such  as  a  current  meter  and  the  flux- 
gate  magnetometer  (Section  1.1).  The  main  purpose  of  these 
filters  is  to  (1)  specify  the  ideal  spectra  densities  in  the 
corresponding  Eq.  (68),  which  are  obtained  from  unit  impulse 
function  inputs  to  the  filters  and  the  use  of  Eq.  (6)  with  T  =  1; 
and  (2)  determine  the  estimators  in  Eq.  (68)  by  using  random 
stationary  series  inputs  to  the  filters  to  generate  the  data 
record  Xg^(t),  Xg2(t)  and,  therefore,  y(t)  for  the  gradiometer 
as  well  as  x^(t)  and  x2(t)  for  the  subsidiary  data  channels. 

The  series  y(t)  simulates  then  the  superposed  time  series  in 
the  main  data  channel,  i.e.,  Xg^(t)  or  the  desired  gradiometer 

signal  and  x„  (t)  or  the  gradiometer  ’’noise";  whereas  the  sub- 
02 

sidiary  time  series  x^(t)  and  x2(t)  denote,  respectively,  the 
internal  wave  measurements  below  the  sea  surface  and  fluxgate 
magnetometer  data  above  the  sea  surface.  Hence,  x^(t)  and  x2(t) 
simulate,  respectively,  the  internal  wave  contribution  to  the 
magnetic  gradients  and  the  background  "noise"  from  conduction 
as  well  as  magnetization  currents  in  nearby  objects. 


r 


"GRADIOMETER” 


FIGURE  2.  Schematic  presentation  of  a  two-input  "gradiometer" 
and  two  subsidiary  channels  for  an  arbitrary  sim¬ 
ulation  of  the  Physical  Dynamics  experiments. 
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The  input  series  n^(t)  and  /  F~ n2(t)  in  Pig.  2  are  random 
stationary  series  that  simulate  uncorrelated  inputs  to  the  G^(f) 
and  G2(f)  gradiometer  filters,  where  the  /F~  parameter  denotes 
noise  in  the  bandwidth  of  the  G2(f)  channel.  The  times  series 
n^(t)  in  the  main  data  channel  is  used  to  simulate  intrinsic 
noise,  where  the  /F~  parameter  denotes  instrument  white  noise 
across  the  whole  frequency  spectrum.  The  output  time  series  y(t) 
of  the  main  data  channel  simulates  then  a  superposition  of  the 
signal  Xg^(t)  with  the  time  series  Xg  (t)  and  n^(t),  both 
of  which  represent  noise  background  including  intrinsic  instru¬ 
ment  noise.  The  effects  of  the  noise  background  in  the  main  data 
channel  can  therefore  be  simulated  by  varying  the  /F~  and  /F^ 
parameters  to  yield  specified  signal  to  noise  ratios.  The  simu¬ 
lated  time  series  for  the  filter  inputs  are  taken  as  stationary 
series  of  random  numbers  with  lengths  of  2048  discrete  points. 

Just  as  was  the  case  of  the  Physical  Dynamics  experiments. 

As  noted  before,  the  above  time  series  [y(t),  x^(t),  x2(t)] 
fix  the  estimators  for  the  spectral  densities  of  the  augmented 
matrix  and,  therefore,  the  estimated  transfer  functions  and  esti¬ 
mated  partial  coherences.  Hence,  by  choosing  the  characteristics 
of  the  foregoing  filters  (or  specified  ideal  spectra)  and  the 
magnitude  of  the  SF~  and  /F^  noise  parameters  (or  signal  to 
noise  ratios),  it  becomes  possible  to  test  the  algorithm  for  the 

/v 

reconstruction  of  the  gradiometer  signal  (S^  ^)  over  a  rather 
wide  range  of  ideal  spectra  characteristics  and  signal- to-noise 
ratios . 

Figure  2  includes  the  reconstruction  of  the  superposed  time 
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series,  y(t),  which  is  as  described  in  Pig.  lb.  The  reconstruction 
process  utilizes  the  filter  output  data  simulating  the  gradiometer 
and  subsidiary  channels  outputs  [y(t),  x^(t),  x2(t)]  to  obtain, 
as  indicated  above,  the  estimators  of  the  spectral  densities,  the 
estimated  transfer  functions  and  the  estimated  gradiometer  signal. 

A  more  detailed  description  of  the  filter  simulation  of  data  channel 
outputs  and  the  reconstruction  process  is  given  in  the  next  sub¬ 
sections.  oh 


3.1  FILTER  SIMULATION  OF  GRADIOMETER  AND  SUBSIDIARY  CHANNELS 


The  determination  of  the  time  series  for  the  filter  outputs 
indicated  in  Pig.  2  [i.e.,  xgi(t),  xg2(t),  y(t),  x^t),  and 
x2(t)]  requires  specifications  of  the  corresponding  transfer 
functions  for  each  filter.  The  gradiometer  and  subsidiary 
channels  are  idealized  by  assuming  that  they  have  time  independ¬ 
ent  properties.  Furthermore,  for  the  gradiometer  output  y(t) 
it  is  assumed  that  the  channels  G1(f)  and  G2(f)  are  characterized 
by  linear  outputs  as  given  in  Eq .  (15).  This  linear  assumption 
implies  two  conditions  (Ref.  2):  (a)  that  the  output  of  several 

simultaneous  inputs  is  equal  to  the  sum  of  outputs  produced  by 
each  input  alone  and  (b)  that  the  output  of  the  product  of  a 
constant  times  an  input  is  equal  to  the  constant  times  the  out¬ 
put  of  the  input  alone.  The  dynamic  characteristics  of  a  con¬ 
stant  parameter  linear  filter  can  be  described  by  a  weighting 
function,  h(x),  which  is  defined  as  the  output  of  the  system  at 
any  time  from  a  unit  impulse  applied  to  it  a  time  t  before.  The 
definition  of  the  weighting  function  h  (t)  for  the  G, (f )  filter, 

for  example,  establishes  both  the  output  xff  (t)  and  transfer 

6] 

function  H  (f)  as 
gl 


and 
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h  (t)  n, (t-T)dx 
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where  Eq.  (74)  is  the  convolution  integral  and  Eq.  (75)  the 

The  lower  limit  of  integration  in 

0  for  t  <  0  (Ref, 


Fourier  transform  of  hg^-r). 


Eq.  (75)  reflects  the  assumption  that  H„  (t) 

gl 


2). 


The  transfer  function  Hg  (f)  is  a  complex  quantity  that  can 
be  expressed  as  Hgi(f)  -  | Hg^Tf )  |  e”1*^,  ,„g 


where  |Hg  (f)|  is 
the  gain  factor  and  $(f)  the  phase  angle.  The  transfer  function 
Hg^for  the  G^(t)  channel  is  given  by  the  Fourier  transform  of 
Eq.  (74),  which  yields 


X  (f)  -  H  (f)  N,  (f )  (76) 

S1  gl  1 

where  N1(f)  =  1  for  a  normalized  random  input  n^t).  A  simula¬ 
tion  of  the  frequency  characteristics  of  the  gradiometer  and 
subsidiary  channels  can  therefore  be  made  by  specifying  the 
transfer  function  for  each  filter,  which  operates  on  the  filter 
input  to  modify  it  according  to  its  frequency  characteristics. 

The  transfer  function  for  the  G^f)  filter  of  the  gradiometer 
would  ideally  simulate  the  desired  spectrum  of  magnetic  field 
gradients  from  internal  waves,  whereas  the  transfer  function 
for  the  K1(f)  filter  would  likewise  simulate  the  spectrum  of 
internal  waves.  However,  for  a  preliminary  test  of  the  algor¬ 
ithm  itself,  the  filter  parameters  can  be  varied  arbitrarily 
to  include  several  combinations  of  filter  characteristics.  The 
problem  of  interest  is  then  to  determine  the  transfer  functions 
for  each  filter,  which  will  yield  specified  ideal  spectra  S  , 

V  V 

S6lgl*  S&2&2’  Sll*  and  S22  for  the  gradiometer  and  subsi¬ 
diary  channels;  i.e.,  with  S^2  =  =  0  in  the  augmented  matrix 

Eq.  (68)  for  the  two  subsidiary  channel  system  because  of  the 
use  of  uncorrelated  inputs.  These  ideal  spectra  can  be  chosen 
by  assuming  their  normalized  shapes,  bandwidths  B,  and  frequencies 
corresponding  to  their  peak  amplitudes,  which,  in  turn,  determine 
the  ideal  transfer  functions  and  H2  [Eq.  (72)]  through  the 
ideal  augmented  matrix  (68). 


The  filter  characteristics  necessary  to  yield  the  specified 
ideal  spectra  can  be  determined  by  using  a  simple,  standard  type, 
digital  recursive  filter  of  the  form  (Ref.  2) 

L 

Cxgl(t)1u=  Cl[nl(t)]u  +Eh£[xgi(t)]u-£  ’  (77> 

1=1 


which  uses  one  input  [n^(t)]u  and  L  previous  outputs  [xg  (t) ]u-£. 
The  summation  term  in  Eq.  (77)  represents  a  discrete,  finite 
sum  equivalent  to  Eq.  (71*)  for  t  =  lAt ,  l  =  1,2  .  .  ,L.  The 
coefficient  h^  represents  L  filter  weights  with  h^  =  h_^  and 
the  subscript  u  -  l  denotes  future  values  of  the  input. 


Equation  (77)  is  illustrated  in  Fig.  3  for  a  second-order 
(L  =  2)  recursive  filter,  where  the  triangles  represent  multi¬ 
plication  by  the  values  within  the  triangles;  the  rectangles 
represent  a  delay  of  At  from  one  point  to  another,  and  the 
circles  represent  summing  operations.  The  Fourier  transform 
of  Eq .  (77 )  yields 


-i2TTf£At 


(78) 


1=1 


where  the  summation  involves  a  polynominal  in  powers  of  exp 
C-i2TrfAt].  From  Eqs.  (76)  and  (78),  the  transfer  function  for 
the  filter  is  given  by 


(f) 


1 


C1 


.  -l2irmt 
n  le 


(79) 


The  determination  of  the  properties  of  the  transfer  function 
Hg  (f)  *s  reduced  to  studying  the  location  and  nature  of  the 
poles  in  the  denominator  of  Eq.  (79),  since  a  filter  is  defined 
to  be  stable  if  all  of  the  roots  of  the  denominator  have  posi¬ 
tive  real  parts. 
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FIGURE  3.  Schematic  representation  of  a  second-order  recursive 
filter  for  gradiometer  and  subsidiary  current  meter 
channels. 


The  application  of  Eq.  (77)  to  each  of  the  four  filters  in 
Fig.  2  yields  the  simulation  of  the  gradiometer  and  subsidiary 
channel  outputs  as  follows: 


\[n‘(t,]» 4  *  c«iK(ll]u.!  (s°> 

A4vt)]u  +  +  c4vt,L2 (si) 

xh  (tl +  xh  (t)lu-i +  xh  (t)]u-2  (8j> 
s[n2  (t)l+  xh  -L  +  st2  ct)lu-2  (83) 


where  the  above  coefficients  may  be  generalized  as  A^,  and 
with  i  *  g^,  g2>  k~2'  These  coefficients,  corresponding  to 
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the  cli,  hli  and  h2i  coefficients  in  a  generalized  Eq.  (11),  can 
be  expressed  in  terms  of  parameters  of  the  specified  ideal  auto- 
spectral  densities.  Using  Eq.  (76),  the  corresponding  ideal  auto- 
spectral  density  functions  S^,  were  arbitrarily  chosen  as 


S1±(f)  -  x±  X»1  =  | H1 | 2 


rcos(7rf/Bi) 
1  +  - S~ - 


(84  ) 


where 


=  008(2*1^/10 

cos  (ttBj, /N) 

(85) 

=  sin(2TrL1/N) 

sin(7rB1/N) 

(86) 

£Li 

N 


(87) 


The  parameter  fQ^  is  the  center  frequency  for  peak  amplitude  of 
the  ith  spectrum,  the  bandwidth  of  the  ith  filter  and  the 
corresponding  center  frequency  parameter  as  given  by  Eq.  (87). 

The  characteristics  of  each  of  the  four  filters  may  then  be 
designated  by  the  choice  of  the  and  parameters.  The  co¬ 
efficients  ,  B^  and  were  then  obtained  by  equating  the 
assumed  form  for  | | 2  in  Eq.  (84)  to  the  corresponding  expres¬ 
sion  obtained  from  a  generalized  Eq.  (79)  for  |H^|2  with  L  =  2 
and  At  =  1/2B.  This  procedure  requires  (1)  the  calculation  of 
XiXi*  in  Eq.  (79),  (2)  equating  the  corresponding  trigonometric 
terms  in  the  denominators  of  Eqs .  (79)  and  (84)  after  putting 
Eq.  (79)  in  a  form  equivalent  to  that  of  Eq.  (84)  by  (a)  retain¬ 
ing  in  Eq.  (79)  only  the  real  parts  and  (b)  dividing  its  result¬ 
ing  numerator  and  denominator  by  c^,  (3)  combining  the  resulting 
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three  equations  so  as  to  solve  c^,  h^  and  h2i  in  Eq.  (79)  in 
terms  of  the  s^  and  r^  variables  in  Eq.  (84)  as  given  by  Eqs. 
(85),  (86);  and  (4)  solving  the  resulting  quartic  equation  in 
terms  of  the  variable  defining  c^,  i.e.,  =  c^/28^.  These 

three  equations  yield  the  A1,  B1  and  Ci  coefficients  in  Eqs. 
(80  to  (83)  as  follows: 


=  c 


li 


=  s 


B,  =  h 


li 


22 
i"i 


i 

4r<  * 


C,  =  h 


2i 


Z±+2 

_  i  g 

2  *1 


where  3^  is  a  root  of  the  quartic  equation 

2 

2i4  *  agi3  “  0gi  -  0*!  +  16  =  0 

where 

a  =  8(s*  +  r*  ) 

8  =  8[4(s’  -  rp  +  1] 

0  =  32  (s*  +  r*). 


(88) 


(89) 


a  root  that  must  also  satisfy  the  condition  of  filter  stability 
given  by 

"hll  -  ^ll  +  11  h2i 

a  h21 

a  condition  that  requires  that  the  root  in  Eq.  (88)  must  be 
real,  positive  and  satisfy 


<  2 


¥■ 


1  +  4  r* 
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Assuming  for  simplicity  that  every  filter  in  Fig.  2  has  the 

same  bandwidth  (B^  =  B)  as  well  as  =  Lg^  +  B  and  -  Lg^  -  B, 

the  ideal  spectra  of  the  gradiometer  and  subsidiary  filters  can 

be  specified  as  illustrated  in  Fig.  4;  i.e.,  from  the  choice  of 

the  filter  parameters  L„  ,  Lff  and  B/2 .  These  parameters  specify 

&1  °2 

also  the  transfer  functions  for  each  filter  (Eq.  79),  and  the 
ideal  transfer  functions,  1^,  H 2  (Fig.  2)  which  can  be  obtained 
either  from  Eq.  (25)  or  (72)  with  S^2  =  S21  =  0,  as 

H,  -  fl£  (90  ) 

S11 

H„  =  ^  .  (91) 

S22 

Hence,  for  a  given  selection  of  filter  parameters,  the  pro¬ 
cedure  used  to  determine  the  ideal  auto-spectra  (sglgl»  Sg2g2, 

Syy,  S.^  and  S22)  as  well  as  the  ideal  transfer  functions  (H^, 

H2)  is  to  (a)  use  the  filter  outputs  to  a  unit-impulse  input  to 
each  filter  so  as  to  generate  the  corresponding  ideal  time  series 
without  variance,  and  (b)  perform  FFT  operations  (e.g.,  Eqs.  5-9) 
with  T  (or  N)  *  1.  For  the  spectral  estimators  as  well  as  esti¬ 
mated  transfer  functions  and  coherences,  the  procedure  consists 
in  (1)  using  the  random  series  inputs  to  the  filters  (Fig.  2) 
so  as  to  obtain  the  output  time  series  (Eqs.  80-83)  as  well  as 
the  gradiometer  output  y(t)  from  Xg^(t)  and  xg2(t)  as  a  function 
of  both  the  filter  parameters  and  signal-to-noise  ratios,  (2) 
performing  FFT  operations  (Eqs.  5-9)  on  the  filter  outputs  for 
y(t),  x^(t)  and  x2(t),  (3)  doing  the  frequency  averaging  of 
each  periodogram  (Eq.  11),  and  (4)  using  the  algorithm  with  the 
smoothed  periodograms  as  spectral  estimators  instead  of  ideal 
spectra. 
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FIGURE  4.  Specified  ideal  spectra  characteristics  for  filter  imple¬ 
mentation  of  gradiometer  and  subsidiary  current  meter 
channels.  The  spectra  Gy,  Ky  are  for  the  signal  channels, 
whereas  the  spectra  G2,  K2  correspond  to  the  noise  chan¬ 
nels. 


3.2  ALGORITHM  FOR  A  SIMULATED  CHANNEL  LINEAR  SYSTEM 

The  basic  results  to  be  obtained  from  the  algorithm  for 
the  reconstruction  of  a  simulated  two-input  and  one-output 
gradiometer  linear  system  (Fig.  2)  are  as  follows:  (1)  the 

A 

estimated*  signal  spectrum  as  given  by  Eq.  (33),  corre¬ 

sponding  to  the  output  Xjr.  (t)  from  the  channel  G, (f),  (2)  esti- 
mated  gain  factors  | H1 |  and  | |  of  the  transfer  functions  for 
the  reconstructed  time  series  that  are  superposed  in  the  main 
data  channel  (Fig.  1),  (3)  the  estimator  for  the  output  of  the 

A  A 

gradiometer,  Syy,  and  (4)  the  estimated  partial  coherence, 
for  the  subsidiary  signal  channel  K1(f)  after  subtracting  in 


*The  estimated  transfer  function  and  partial  coherences  are  not 
estimators,  because  they  are  based  only  on  the  use  of  the  esti¬ 
mators  for  their  spectral  densities  and  not  on  the  averaging 
of  themselves.  For  example,  <H^>  /  >/<§^>,  but  lHil  =  <|H^|> 
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the  main  reconstructed  gradiometer  channel  the  predictions  from 
the  ^(f)  channel.  Each  of  the  foregoing  results  is  given  by 
the  algorithm  in  terms  of  the  estimators  for  the  appropriate 
auto-  and  cross-spectral  density  functions  in  the  augmented 
matrix  given  in  Eq.  (68).  However,  in  order  to  facilitate  the 
Interpretation  of  results  from  subsequent  numerical  simulations, 
it  is  of  interest  to  bring  out  the  relationships  among  the  esti¬ 
mators  and  estimated  .parameters  directly  from  the  subsequent 
brief  considerations. 


3.2.1.  Estimated  Gradiometer  Signal  Spectrum 


As  indicated  in  Fig.  2,  the  estimated  gradiometer  signal 
spectrum  is  conceptually  obtained  from  the  time  series  5^(t)j 
which  is  an  output  reconstructed  by  the  transfer  function  de¬ 
scribing  the  contribution  of  the  time  series  x^(t)  or  internal 
waves  to  the  magnetic  gradients.  Hence,  Eqs.  (74)  and  (75)  be¬ 
come  applicable  to  the  transfer  function  between  C1(t)  and  x^t); 
i.e.,  by  replacing  xg^(t),  hg^x)  and  n^t)  by  ^(t),  h-^x)  and 
x1(t),  respectively.  For  a  pair  of  times  t  and  t  +  x,  the  pro¬ 
duct  ^(t)  ^(t+x)  is  given  from  Eq.  (74)  by 


Cx(t) 


5x(t+x) 


o 

7/ 


h(Oh(n) 


x1(t— ^ ) 


x^  (t+x— ri)d£dn 


(92) 


Since  the  auto-correlation  function  for  a  stationary  series  is 
given  by 

/'T 

dt  (93  ) 


R 


-  ,  (x)  =  lim  1  f  ^(t)  £  (t+x) 

T-“  T  J  ±  ± 


the  combination  of  Eqs.  (92)  and  (93) >  together  with  Eq.  (1), 
yields  the  important  result 


5,  rr  (f)  =  |H.(f)|2  Sy  (f)  =  |  H,  (  f )  | 
*1*1  A  X1X1 


Sll(f) 


(94  ) 
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Equation  (94)  indicates  that  the  estimated  spectrum  of  the 

/v 

signal,  ^  ,  is  Siven  by  the  product  of  the  square  of  the 
estimated  gain  factor  of  the  transfer  function  and  the 
estimator  of  the  spectrum  S.,.;  i.e.,  the  estimated  signal  spec- 
trum,  is  determined  from  the  corresponding  values  of  | | 

and  S.^.  The  gain  factor  |H.J,  for  example,  is  based  on  the 
estimators  for  the  auto-  and  cross-spectral  densities  of  y(t), 
x1(t)  and  x2(t)  as  a  function  of  frequency  averaging  and  signal- 
to-noise  ratios,  since  y(t)  is  specified  by  the  /F  and 
parameters  in  the  inputs  to  the  gradiometer  filters  (Fig.  2). 
Hence,  IH^J  is  based  on  the  full  augmented  matrix  (68)  i.e., 
where  S12  and  S21,  for  example,  are  nonzero  due  to  the  bias  in 
the  estimators  from  the  averaging  of  the  periodograms .  Hence, 
Eq.  (94)  will  be  an  important  criterion  to  test  the  bias  of  the 


estimator  Sg  ^  from  its  comparison  with  the  corresponding  value, 
Sg  ,  which  is  obtained  from  the  G^(f)  filter  output  in  the 
time  domain  to  a  unit  impulse  input  and  the  subsequent  applica¬ 


tion  of  Eq.  (6)  with  T  =  1. 


3.2.2  Estimated  Gradiometer  Transfer  Functions 

Since  the  linear  assumption  given  by  Eq .  (15)  indicates 
/\ 

that  the  output  y(t)  for  the  gradiometer  (Fig.  2)  can  be  ex¬ 
pressed  by  the  sum  of  the  individual  outputs  £^(t),  then 

2 

y(t)  “X>i(t)  (95) 

i-1 


where  £^(t)  is  defined  as  that  part  of  the  output  which  is  pro¬ 
duced  by  the  ith  input  when  all  the  other  inputs  are  zero. 
Expressing  each  output  E1(t)  in  terms  of  its  respective  input 
x^(t)  through  the  convolution  integral  [Eq.  (74)], 

/s  2  r  » 

y(t)  »  ^  /  h1(x)  x1(t-x)dT.  (96) 

i=l  J  o 
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Calculation  of  the  cross-correlation  function,  R  (t),  in  terms 

xy 

of  the  output  as  given  by  Eq.  (95),  in  a  manner  similar  to 
that  leading  to  the  results  given  by  Eqs.  (92)  and  (9*0  for  a 
single  channel  yields 


which  yields 


siy  =  f  Vli 

(  97) 

J  =  1 

ly  =  H1S11  +  H2S12 

(  98) 

2y  =  H1S21  +  H2S22 

(99) 

The  input  auto-  and  cross-spectral  density  functions  in  Eqs.  (58) 
and  (99)  can  be  used  to  define  the  input  coherence 


(100  ) 


where  the  numerator  is  the  absolute  value  of  the  complex  cross- 
spectral  density  function  for  the  two  inputs  x1(t)  and  x2(t) 
and  the  denominator  the  corresponding  real-valued  autocorrela¬ 
tion  density  functions.  The  input  coherence  in  the  frequency 
domain  is  thus  equivalent  to  the  normalized  cross-correlation 
coefficient  in  the  time  domain,  i.e.,  both  are  real-valued 
functions  with  magnitudes  between  zero  and  unity.  Just  as  the 
cross-correlation  coefficient  establishes  the  degree  of  corre¬ 
lation  between  two  given  time  series  in  the  time  domain,  the 
input  coherences  do  likewise  in  the  frequency  domain.  Hence, 
for  the  two  uncorrelated  times  series  x^t)  and  x2(t)  represent¬ 
ing  internal  waves  and  background  noise  from  conduction  and  mag¬ 
netization  currents,  both  the  cross-correlation  function  and 
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input  coherence  must  be  zero.  Eqs.  (98)  to  (100)  yield  the 
transfer  functions  for  the  reconstructed  gradiometer  channels  as 


„  =  Sly  ~  S12S2y/S22  =  Sly  ~  S12S2y/S22 

1  S11  '  S12S21/S22  sii  (1  _ 

and 

h  =  S2y  ~  S21Sly/Sll  =  S2y  ~  S21Sly/Sll 

2  S22  "  S21S12/S11  S22  (1  -  Y22) 


(101) 


(  102) 


Eqs.  (101)  and  (102)  will  be  used  to  obtain  the  estimated  trans¬ 
fer  functions  and  gain  factors  in  terms  of  the  estimators  for 
the  relevant  auto-  and  cross-spectral  densities  functions. 

These  equations  also  yield  the  corresponding  ideal  reference 
values  [Eqs.  (90,  91)]  for  which  Sl2  =  S21  *  y22  =  0. 

3.2.3  Estimator  for  the  Gradiometer  Output  Spectra 

y  /v 

The  estimator  for  the  gradiometer  output  spectrum,  S  ,  is 

«y  v 

obtained  by  a  procedure  similar  to  that  leading  to  Eq.  (97),  ex¬ 
cept  that  Instead  of  Eqs.  (92)  and  (93)  for  the  auto-correlation 
function  of  the  output  E^t),  it  becomes  necessary  to  use  the 
auto-correlation  for  the  output  y(t).  The  result  is  then  given 

by 

syy  ”  £  Z  hJ  Hj  su  ( 103) 

i-1  J-l 

where  represents  the  cross-spectral  density  function  be¬ 
tween  the  inputs  x^(t)  and  Xj(t).  Eq.  (103)  then  yields 

syy  -  lHll2su  +  Hj  H2S12  +  H*  H1S21  +  |H2|2S22  .  (1M) 

A 

Eq.  (104)  shows  the  relationship  among  the  estimators  for  S- 
and  those  involving  the  gain  factors,  transfer  function  and 
spectral  density  functions  derived  from  the  two  current  meters 
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in  Pig.  2.  Although  for  the  ideal  case  S12  =  s2i  =  °>  the  ideal 
value  of  S  was  determined  by  computing  the  two  remaining  terms 

«/  J 

in  Eq.  ( 1 0 4 )  from  an  addition  of  the  ideal  values  of  Sg-^g.^  and 
sg2S2* 

3.2.4  Estimated  Partial  Coherence 

Tne  partial  coherence  is  given  by  a  ratio  similar  to  that 
for  the  input  coherence  [Eq.  (100)],  except  that  the  power  spec¬ 
tral  density  functions  involve  the  residual  Ax^(t)  and  Ay(t) 
time  series  obtained  by  subtraction  from  x, (t)  and  y(t)  their 
respective  linear  least-square  predictions  (x-^y)  from  x2(t), 
i.e.  , 

Axl  =  xi^^  ~  xx(t)  =  x1(t)  -  (Rx  x  /Rx  x  )  x2(t)  (105) 


Ay  *  y(t)  -  y (t )  =  y(t)  -  (R  /R  )  x?(t) 

a K2  2  ^ 


c  106) 


where  R*^,  Rx2Xl  and  Rx2y  denote  the  respective  cross-correlation 
functions  for  the  indicated  inputs  x.^  and  output  y.  The  partial 
coherence  for  the  residuals  Ax-^  and  Ay  time  series  is  then  given 
by 


Ax1Ay 


(f)  *  vjy.2  = 


S11.2  Syy . 2 


(  107) 


where  the  subindex  notation  has  again  been  simplified  by  dropping 
the  respective  x's  and  the  symbol  .2  denotes  subtractions  of  the 
respective  predictions  from  x2(t)  as  indicated  by  Eqs.  (105) 
and  (106).  The  importance  of  the  partial  coherence  comes  about 
from  the  fact  that  its  value  must  be  unity  at  every  frequency 
[Eqs.  (^6),  (60),  and  (6l)],  a  result  that  is  also  obtained  from 
expressions  for  the  numerator  and  each  term  in  the  denominator 
of  Eq.  (107),  as  given  in  the  sequel. 


The  spectral  density  functions  necessary  for  the  denominator 
of  Eq.  (107)  are  obtained  by  a  procedure  similar  to  that  lead¬ 
ing  to  Eq.  (97),  except  that  consideration  must  instead  be  given 
to  the  residual  time  series  given  by  Eqs.  (105)  and  (106).  If 

A 

instead  of  Eq.  (74),  the  least  square  prediction  y(t)  from  x2(t) 
is  given  by 


y(t) 


hQ(T)x2(t-T)dT 


(  108) 


where  h  (t)  is  a  weighting  function  determined  so  that  the  mean 
0  2 

square  error  e  for  the  expected  value  of  the  auto-correlation 
of  Ay  is  a  minimum,  the  result  is  then  a  relationship  between 
the  cross-  and  auto-correlation,  whose  Fourier  transform  yields 
Hq  =  S2y/S22.  Subsequent  calculation  of  the  auto-correlation 
function  RAyAy(T)  and  its  Fourier  transform  leads  to  the  result 


Ay  Ay 


=  S 


yy-2 


'yy 


S2ySy2 

S22 


( 109) 


where  the  subscript  2  denotes  the  series  x2(t)  from  which  the 
prediction  y  is  made.  Equation  (109)  applies  also  to  the  other 
spectral  density  functions  in  the  denominator  of  Eq.  (107),  i.e. 


11.2 


=  S 


Ax^Ax^ 


=  S 


’12’ 
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-  sn(1  - 
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The  numerator  of  Eq.  (107)  is  likewise  obtained  by  a  procedure 
similar  to  that  leading  to  Eq.  (109),  except  that  the  cross- 
correlation  function  R^x  ^(t)  must  be  utilized  instead  of  R^y^y(T)- 
The  result  is  then  1 


S  =  S  =  S  S12S2y 

ly.2  -  bAX;LAy  sly  S22  * 

Substituting  Eqs.  (98)  and  (99)  for  Sly  and  S2y  in  Eq. 
there  follows 


(111) 

(99) 
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(112  ) 


’ly.2 


=  H 


1  11 


(1_Y12 ) 


for  the  numerator  of  Eq.  (107).  Similarly,  substituting  Eqs. 
(99)  for  S2y  and  (104)  for  Syy  into  Eq.  (109)  it  is  obtained 


yy 


2  =  |h1|2s11(i-y^2) 


(113) 


Hence,  substitution  of  Eqs.  (110),  (112)  and  (113)  into  (107) 
yields 

*ly.2  *  1  <U4> 

at  any  frequency,  a  result  that  holds  for  z(t)  ~  0  in  Eq.  (12) 
whether  or  not  the  two  inputs  x^(t)  and  x2(t)  are  correlated. 
Hence,  an  Important  question  in  the  test  of  the  algorithm  is 
the  number  of  frequencies  (M)  required  in  the  averaging  of  the 
periodograms  [e.g.,  Eqs.  (6),  (8)]  so  as  to  reproduce  the  given 

A 

signal,  i.e.,  to  obtain  «  Sg^-^as  well  as  the  behavior  of 

the  input  (y22  0)  and  partial  (y|y  2  =  1)  coherences  at  every 

frequency. 


3.3  SCOPE  OF  NUMERICAL  SIMULATIONS 

Figure  5  shows  the  specified,  normalized  ideal  spectra 
(SgiEi*  Sg2g2>  S22)  *'or  ^our  cases  °f  a  two-input  gradiom- 

eter  and  two  subsidiary  single-input  channels  as  illustrated  in 
Fig.  2.  For  a  given  case,  each  channel  is  identifiable  by  the 
frequency  parameter  for  peak  amplitude  of  the  ideal  spectra  as 
follows  (Fig.  4):  (a)  Sg^g^  at  Lgl ,  with  a  bandwidth  B  at  the 

half-power  level,  (b)  Sg2g2  at  Lg2,  (c)  at  Lgl  +  B  and  (d) 
S22  at  Lg2  -  B.  The  ideal  spectra  for  the  gradiometer  output, 

S  ,  is  given  by  the  sum  of  the  corresponding  spectra  for  each 

y  y 

channel  (sglgl,  Sg2g2)  of  gradi°rceter .  The  frequency  range 
corresponds  to  one-half  the  length  of  the  discrete  stationary 
series  input  consisting  of  2048  points  in  the  time  domain.  The 
four  cases  shown  in  Fig.  5  are  characterized  by  both  the  center 
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FIGURE  5.  Normalized  amplitudes  of  specified 
Ideal  spectra  for  gradlometer  and 
subsidiary  channels  as  a  function 
of  frequency. 
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frequency  parameter  for  peak  amplitude  and  the  bandwidth  of  the 
spectra  (Pig.  4).  The  first  two  cases  are  for  low  center  fre¬ 
quency  parameters  and  narrow  bandwidths,  whereas  the  last  two 
consider  broad  bandwidths.  Furthermore,  the  first  low-frequency 
case  is  characterized  by  a  ratio  Lg2/Lg^  >  1,  whereas  the  second 
one  has  Lg2/Lg^  —  1.  The  two  cases  for  high  frequency  have 
Lg2/Lgi  ^  1  but  somewhat  different  spectra. 

Table  1  gives  the  values  of  the  center  frequency  parameters 
for  peak  amplitude  of  the  spectra  (Lg  ,  Lg2)  and  half  the  band¬ 
width  (B/2)  for  each  of  the  four  cases  in  Fig.  5.  It  should  be 
noted  that  half  the  bandwidth  is  the  reference  for  the  number 
M  of  frequencies  to  be  used  at  either  side  of  a  given  frequency 
in  the  smoothing  out  of  the  periodograms  for  each  spectral 
density.  As  indicated  in  Eq.  (11),  the  degree  of  frequency 
averaging  in  a  given  periodogram  is  given  by  2M  +  1  or,  approxi¬ 
mately,  by  M/(B/2)  relative  to  the  bandwidth  of  the  filters; 
i.e.,  a  value  of  M/(B/2)=1  indicates  frequency  averaging  over 
a  window  equal  to  the  bandwidth  of  the  specified  spectra.  Table 
1  includes  the  range  of  the  gradiometer  signal-to-noise  ratio 
(S/Np)  and  the  corresponding  values  of  the  /F~  parameter  shown 
in  Fig.  2.  The  gradiometer  signal-to-noise  ratio  is  inversely 
proportional  to  F  and  it  was  varied  in  the  range  1/10  to  10, 
which  yielded  the  corresponding  /F~  value,  shown  in  Table  1. 

The  table  includes  also  the  range  of  the  signal-to-instrument 
noise  ratio  (S/N^  and  the  corresponding  values  of  the 
parameters  indicated  in  Fig.  2.  The  white  noise  power  within  the 
bandwidth  of  the  filter  is  given  by  the  corresponding  fraction 
of  the  white  noise  across  the  whole  frequency  spectrum.  There¬ 
fore,  the  signal-to-instrument  noise  ratio  is  inversely  propor¬ 
tional  to  the  product  of  F^  and  the  fraction  of  the  area  within 
the  bandwidth  signal,  i.e.. 


FT  (A  /102K7 

1  8lsl 


(115) 
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where  Ag  g  denotes  the  area  under  the  ideal  spectrum  Sg^g^  and 
1024  is  the  frequency  range  (T/2)  for  the  input  time  series  of 
2048  points.  The  signal-to-instrument  noise  ratio  values  were 
taken  as  °°,  10  and  1;  which  yielded  the  corresponding  /F^ 
values  in  Table  1.  These  values  of  the  signal-to-noise  instru¬ 
ment  ratio  represent,  respectively,  (a)  ideal  conditions  given 
by  z(t)  =  0  in  Eq.  (12),  (b)  signals  with  a  moderate  instrument 
noise  background,  and  (c)  signals  that  are  expected  to  be  with¬ 
in  the  threshold  of  detectability. 


TABLE 

1  .  FILTER  PARAMETERS 

AND  SIGNAL-T0- 

NOISE  RATIOS 

Case 

L9l 

l2i 

B/2 

s/nf 

77 

S/N. 

I 

50 

220 

10 

1 

1 

00 ,  10, 

1 

0,  1.84,  5.82 

II 

50 

60 

20 

1 

1 

“*  10, 

1 

0,  1.43,  4.51 

III 

500 

500 

100 

1 

1 

®,  10, 

1 

0,  0.59,  1.87 

IV-a 

500 

600 

100 

1 

1 

°°»  10, 

1 

0,  0.58,  1.85 

IV-b 

500 

600 

100 

10 

0.316 

OO  -- 

-- 

0,  -  ,  - 

IV-c 

500 

600 

100 

1/10 

3.162 

— , 

— 

0,  -  ,  - 

The  filter  parameters  in  Table  1  define  the  ideal  transfer 
functions  of  the  simulated  gradiometer  and  subsidiary  channels 
in  Fig.  2.  The  use  of  uncorrelated  input  time  series  [n1(t), 
n2(t)]  generate  the  output  time  series  given  by  Eqs. (80)— (83) 
as  well  as  the  output  time  series  y(t),  which  includes  /F^n^ 
(Fig.  2).  The  subsequent  use  of  equations  such  as  Eqs.  ( 5 ) — ( 9 ) 
gives  the  corresponding  periodograms ,  which  are  then  averaged 
with  respect  to  frequency  as  illustrated  by  Eq.  (11).  The  es¬ 
timators  for  the  spectral  density  functions  (i.e.,  smoothed 
periodograms)  yield  then  the  estimated  transfer  functions  [Eqs. 
(101),  (102).],  the  reconstructed  spectra  for  the  gradiometer 
output  [Eq.  (104)J,  the  imbedded  signal  in  the  gradiometer  [Eq. 
(9-4)J,  and  the  corresponding  partial  coherence  [Eq.  (107),  (109- 
(111)]. 
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3.4  INTERPRETATION  OF  NUMERICAL  SIMULATIONS 


Figure  2  indicates  that  the  two  subsidiary  channels  do 
not  process  stationary  time  series  with  instrument  noise.  Hence, 

/\  /v 

the  estimators  for  the  spectra  and  S^2  as  well  as  the  esti¬ 
mated  input  coherence  y ^  2  will  be  shown  along  with  other  data 
only  for  one  signal-to-instrument  noise  ratio. 

Figure  2  also  shows  that  the  simulated  gradiometer  output 
time  series  xg^(t)  and  Xg2(t)  are  unaffected  by  the  instrument 
noise  in  the  main  data  channel.  Hence,  only  auto-  and  cross¬ 
spectra  densities  involving  the  main  data  channel  output,  y(t), 
will  be  modified  by  the  /Fj  white  noise  parameter.  Likewise,  the 
output  time  series  Xg-^t)  is  unaffected  by  the  gradiometer  noise 
parameter  /F~ .  Hence,  only  the  auto-spectral  density  as 

well  as  the  auto-  and  cross-spectral  densities  involving  the 
main  data  channel  output  y(t)  with  or  without  instrument  noise 
will  be  modified  by  the  /F  parameter.  It  is  important  to  empha¬ 
size  that  the  same  set  of  random  stationary  series  was  always 
utilized  for  n^Ct),  n2(t)  and  n^(t),  which  were  characterized 
by  cross-correlation  coefficients  that  were  very  small.  Further¬ 
more,  each  time  series  was  assumed  to  be  given  by  a  sequence  of 
random  numbers  of  2048  points,  and  their  frequencies  to  be  rep¬ 
resented  by  the  order  of  each  term  in  such  series  in  the  fre¬ 
quency  domain,  i.e.,  by  1,  2,  ....1024. 

Before  presenting  the  data  for  each  case,  it  is  important 
to  indicate  some  typical  trends  in  the  periodograms .  For  rela¬ 
tively  unaveraged  conditions  (say  M  =  1),  the  input  coherence 
~2 

Y  is  highly  noisy,  i.e.,  it  varies  as  white  noise  between  zero 
and  unity.  These  large  fluctuations  are  smoothed  out  as  the 
number  of  frequencies,  M,  in  the  frequency  averaging  (2M  +  1) 
increases  toward  the  bandwidth  value,  B,  of  the  specified  spectra 
in  Fig.  5;  i.e.,  as  M/(B/2)  -*■  1.  Similar  trends  take  place  for 
the  gain  factors,  IH^J,  and  some  spectral  densities  when  the 
signal-to-instrument  noise  is  low,  i.e.,  their  magnitudes  can 
become  extremely  high  for  relatively  unaveraged  (M  =  1)  periodo¬ 
grams  . 
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4.  RESULTS  OF  NUMERICAL  SIMULATIONS 


The  numerical  simulation  of  a  two-input  gradiometer  and  two 

subsidiary  single-input  channels  (Pig.  2)  considered  the  ideal 

spectra  shown  in  Fig.  5»  i.e.,  (a)  two  cases  with  low  values  of 

the  center  frequency  parameters  L  and  L_  (Eq.  87)  and  (b)  two 

S1  62 

cases  with  high  values  of  the  center  frequency  parameters  together 
with  broader  bandwidths.  The  low  frequency  cases  are  further¬ 
more  characterized  by  two  values  of  the  Lg  /Lg  ratio,  i.e., 
when  the  gradiometer  "noise"  channel  has  a2cen^er  frequency  for 
peak  amplitude  of  the  spectrum  that  is  significantly  higher 
than  that  for  the  gradiometer  signal  channel  and  when  these 
center  frequencies  are  nearly  equal. 

As  indicated  in  the  preceding  section,  the  results  of  the 
numerical  simulations  for  each  case  put  emphasis  on  (a)  the 

/v 

estimated  gradiometer  signal  spectrum  (Eq.  9*0,  (b)  the 

estimated  gain  factors  |H^|  and  |H2|  of  the  simulated  gradiom¬ 
eter  transfer  functions;  the  former  appears  in  the  expression 

A 

for  while  both  do  so  in  that  for  the  estimator  of  the 

gradiometer  output  spectrum  S  (Eq.  104),  (c)  the  estimators 

yy  ~  ~ 

for  each  spectrum  of  the  subsidiary  channels,  and  §22’’ 

again,  the  former  appears  in  the  expression  for  Sp1 f, ,  whereas 
both  do  so  in  that  for  S  ,  (d)  the  estimator  S  ,  (e)  the 

v  j  <y  v 

estimated  partial  coherences  for  the  gradiometer  signal  channel 
(Eq.  107),  and  (f)  the  estimated  input  coherence  (Eq.  100). 

These  results  are,  in  general,  given  for  each  case  as  a  function 
of  both  the  number  of  frequencies  (M)  used  in  the  averaging  of 
the  periodograms  (Eq.  11)  and  each  of  three  characteristic  sig- 
nal-to-instrument  noise  ratios  (Table  1).  As  noted  previously, 
the  ratio  M/(B/2)  indicates  the  length  of  the  averaging  window 
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relative  to  the  bandwidth  of  the  specified  ideal  spectra.  The 
three  characteristic  signal-to-instrument  noise  ratios  corre¬ 
spond  to  an  ideal  case  of  no  instrument  noise  (S/N^  =  “),  de¬ 
tectable  signal-to-noise  ratios  (S/N^  =  10)  and  signals  that 
are  within  the  threshold  of  detectability  (S/N^  =  1). 

The  results  for  each  foregoing  case  and  parameter  include 
comparisons  of  the  ideal  and  estimated  or  estimators  parameters, 
so  as  to  exhibit  the  bias  in  the  estimated  or  estimator  param¬ 
eters  as  a  function  of  the  frequency  averaging,  M/(B/2),  for  the 
three  typical  signal-to-noise  ratios.  The  ideal  values  corre¬ 
spond  to  the  filter  outputs  (Fig.  2)  to  unit  impulse  inputs, 
whereas  the  estimated  and  estimators  parameters,  respectively,  are 
for  the  reconstructed  outputs  of  the  gradiometer  and  actual  out¬ 
puts  of  the  subsidiary  channels  (Figs,  lb  and  2).  The  results 
are  then  plotted  as  a  function  of  simulated  frequency,  i.e.,  by  the 
order  of  the  random  stationary  terms  in  the  input  time  series 
to  the  filters  (Fig.  2). 

4.1  LOW-FREQUENCY  GRADIOMETER  WITH  lQ  /La  >  1 

y2 

Figures  6-27,  in  the  subsequent  subsection,  describe  the 
results  for  the  foregoing  parameters  for  Case  I  in  Fig.  5,  which 
is  characterized  by  L  =  50,  L  =  220,  L  /L  *  4.11,  and 

,  6i  6?  g2  B1 

B/2  -  10.  1  d  d  1 


4.1.1  Estimated  Gradiometer  Signal  Spectrum 

Figures  6-8  show  comparisons  of  ideal  and  estimated  spectra 
for  the  gradiometer  signal  channel  as  a  function  of  frequency 
averaging  and  signal-to-instrument  noise  ratio  for  S/N.^  =  “,  10, 
and  1,  respectively.  These  figures  indicate  the  following: 

A 

•  Figure  6  compares  the  ideal  (S,,  e  )  and  estimated  (Sc  c  ) 

6161  ^1*1 

signal  spectra  for  ideal  conditions  of  no  instrument 
noise  and  K/(B/2)  that  are  as  high  as  twice  the  bandwidth 
of  the  specified  spectrum.  Although  not  shown  in  these 
figures,  the  relatively  unaveraged  (M=l)  noisy  values  for 

A 

Sr  f  show  a  peak  amplitude  that  is  larger  than  the 
*1*1 


ideal  peak  value  at  Lg^  =  50  by  a  factor  of  about  1.45, 
This  figure  then  indicates  that  frequency  averaging 
smooths  out  the  fluctuations  about  the  mean  signal,  but 
it  underestimates  significantly  the  actual  spectrum  even 
when  M  =  5  or  M/(B/2)  =  0.5,  i.e.,  when  the  length  of 
the  averaging  window  equals  about  half  the  bandwidth  of 
the  specified  signal. 

A 

•  Figure  7  compares  the  ideal  (Sglgl)  and  estimated  (S^^) 
signal  spectra  when  the  signal-to-instrument  noise  ratio 
S/N.  =  10.  Because  of  the  increased  noise  relative  to 
the  previous  figure,  the  estimated  spectrum  tends  toward 
the  ideal  value  when  1  <  M/(B/2)  <  2.  Although  not  shown 

A 

in  these  figures,  the  values  of  S^^  for  M  =  1  are  highly 
noisy  at  every  frequency,  and  the  magnitude  of  S^^  in 
the  neighborhood  of  Lg^  =  50  is  as  high  as  10  times  that 
of  the  ideal  Sglgl  values.  Figure  7  for  detectable 
signal-to-instrument  noise  ratios  indicates  that  an  un¬ 
biased  estimate  of  the  peak  amplitude  of  the  signal 
Sg^i  wou-1-d  require  an  averaging  window  given  by  1  < 
M/(B/2)  <  2.  However,  the  frequency  characteristics  of 
the  estimated  signal  would  not  be  clearly  detectable  due 
to  the  distortions  of  the  shape  of  the  estimated  signal 
by  the  smoothing  out  of  the  periodogram. 

A 

•  Figure  8  compares  the  ideal  (Sg^i)  and  estimated  (S^i^) 
signal  spectra  for  S/N^  =  1,  i.e.,  when  the  signal-to- 
instrument  noise  ratio  is  within  the  threshold  of  detect¬ 
ability.  This  figure  indicates  even  greater  overesti- 

A 

mates  in  the  peak  amplitude  of  S^^  at  Lgi  =  50  for 
M/(B/2)  =0.5  than  those  in  the  previous  figures  for 
S/N.  =  «  and  10.  Although  not  shown  in  this  figure, 
the  magnitude  of  S^^  for  M  =  1  in  the  neighborhood  of 
Lgi  =  50  is  as  high  as  100  times  that  of  the  ideal  Sg^i 
value.  The  results  in  Fig.  8  for  S/Ni  =  1  indicate  that 
an  unbiased  estimate  of  the  peak  amplitude  of  the  signal 
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FIGURE  6.  Comparison  of  ideal  and  reconstructed  signal  spectr 
for  low-frequency  filter  with  =  4.4  for  S/N 
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FIGURE  7.  Comparison  of  Ideal  and  reconstructed  signal  spectra 
for  low-frequency  filters  with  Ugg/Lg-j  *  4.4  for 
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Comparison  of  ideal  and  reconstructed  signal  spectra 
for  low-frequency  filters  with  Lg9/Lgi  *  4.4  for 


Sg  g  would  require  frequency  averaging  given  by  M/(B/2) 

Z  2.  Furthermore,  the  detectability  of  even  the  peak 
amplitude  of  an  unknown  Sg^g^  spectrum  would  be  difficult 
due  to  the  comparable  white  noise  at  other  frequencies. 

4.1.2  Estimated  Gain  Factors 

Figures  9-13  and  14-18  show  the  estimated  gain  factors  for 
the  gradiometer  signal  and  "noise'’  channels,  respectively,  as  a 
function  of  frequency  averaging  and  the  signal-to-instrument 
noise  ratio  S/Ni  *  »,  10,  and  1.  These  figures  indicate  the 
following : 

A 

•  Figure  9  for  the  estimated  gain  factor  Ih^I  for  the 
gradiometer  signal  channel  and  S/J^  =  ®  indicates  that 
the  transfer  function  itself  would  be  detectable  even 
without  averaging,  i.e.,  for  average  values  as  low  as 
M/(B/2)  =  1/10  or  0.1.  However,  due  to  numerical  noise 
from  the  algorithm,  a  spurious  signal  appears  in  the 
neighborhood  of  Lg2-  It  is  to  be  noted  that  increasing 
M/(B/2 )  as  required  to  detect  Sglgl  not  only  introduces 
a  corresponding  bias  in  the  amplitude  of  the  gain  factor 
but  also  magnifies  the  spurious  signal  at  Lg2.  The 
reason  for  this  spurious  signal  is  the  fact  that  in  the 
neighborhood  of  this  frequency,  the  denominator  in  Eq. 
(101)  decreases  faster  than  the  numerator. 

A 

,  Figure  10  for  |H1|  and  S/Nj  =  10  indicates  drastic  in¬ 
creases  in  the  amplitude  of  the  gain  factor  at  fre¬ 
quencies  higher  than  Lg2-  Although  not  shown  in  this 
figure,  the  amplitude  of  |H^|  for  M  ■  1  is  extremely 
noisy  at  frequencies  f  >  400,  as  indicated  by  values 

A 

1^1  >  1000.  The  significance  of  the  high  values  of 
the  gain  factor  at  the  high  frequencies  is  that  they 

A  A 

contribute  to  the  white  noise  of  (Fig.  7)  and  Syy, 

as  will  be  evident  later.  Figure  11  compares  the  ideal 
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and  estimated  gain  factors  as  a  function  of  M,  which 
indicates  that  M/(B/2)  *  2  in  order  to  obtain  unbiased 
estimated  gain  factors  at  Lg^.  However,  the  numerical 
noise  at  frequencies  higher  than  Lg^  would  make  the  gain 
factor  for  the  signal  gradiometer  channel  undetectable. 

A 

•  Figure  12  for  | |  and  S/I^  =  1  indicates  larger  ampli¬ 
tudes  for  the  gain  factor  by  a  factor  of  about  three 
than  those  for  S/Ni  =  10  in  Fig.  10.  Figure  13  compares 
the  ideal  and  estimated  gain  factors  as  a  function  of  M, 
which  shows  trends  somewhat  similar  to  those  in  Fig.  11 
for  S/N.^  =  10. 

A 

•  Figure  14  for  the  estimated  gain  factor  |H2|  for  the 
gradiometer  "noise"  channel  and  S/N,  =  »  shows  trends 

a  X 

similar  to  those  for  | |  in  Fig.  9  for  the  same  signal- 
tc-instrument  noise  ratio,  except  that  the  numerical 

noise  takes  place  now  at  the  L-  frequency. 

*  ^1 

•  Figure  15  for  | |  and  S/Ni  «  10  shows  trends  similar  to 
those  for  | H1 |  in  Fig.  10  for  the  same  signal- to- instru¬ 
ment  noise  ratio,  except  that  the  amplitude  of  the  |H2| 
are  lower  by  a  factor  of  about  2.7  for  M  =  5,  a  factor 
that  becomes  somewhat  larger  as  M  increases.  Figure  16 
compares  the  ideal  and  estimated  gain  factors  as  a  func¬ 
tion  of  M,  and  shows  trends  similar  to  those  in  Fig.  11 

/v 

for  | H2 |  at  the  same  signal-to-instrument  noise  ratio. 

A 

•  Figure  17  for  | H0 I  and  S/N^  =  1  indicates  a  trend  similar 
to  those  in  Fig.  12  for  |H1|  and  same  signal-to-instru¬ 
ment  noise  ratio,  except  that  the  amplitudes  of  |H0|  are 

A  U 

lower  than  those  for  |H, I  by  a  factor  of  about  1/2. 

Figure  18  shows  comparisons  of  ideal  and  estimated  |H2| 
gain  factors  for  S/N.  =  1,  and  shows  trends  similar  to 
those  in  Fig.  13  for  |  |  and  same  signal-to-instrument 

noise  ratio. 
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FIGURE  9, 
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Comparison  of  Ideal  and  estimated  gain  factors 
for  gradlometer  signal  channel  for  low-frequency 
filters  with  Lg^/Lg^  =  4.4  for  S/N^  =  ®. 
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FIGURE  11.  Comparison  of  Ideal  and  estimated  gain 
for  gradlometer  signal  channel  for  low 
filters  with  Lg2/Lgj  *  4.4  for  S/N1  «= 
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FIGURE  12.  Estimated  gain  factor  for  gradlometer  signal 
channel  for  low-frequency  filters  with 
L92/L9l  s  4-4  for  S/Ni  *  !• 
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CASE  I.  S/Mf  - 1,  S/N)  =  1 


UM7MI 


FIGURE  13.  Comparison  of  Ideal  and  estimated  gain  factors  for 
gradiometer  noise  channel  for  low-frequency  filters 
with  Lg 2 / Lg i  s  4.4  for  S/Ni  =  1. 
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4.1.3  Estimators  for  Grad lometer  Output  Spectra 

Figures  19-23  show  comparisons  of  the  ideal  and  estimated 
spectra  for  the  subsidiary  channels  and  gradiometer  output 
(Fig.  2)  as  a  function  of  frequency  averaging  and  signal-to-in- 
strument  noise  ratio  for  the  gradiometer  output.  These  figures 
indicate  the  following: 

A 

•  Figure  19  for  and  any  S/N^  value  shows  that  signifi¬ 
cant  bias  is  introduced  for  averaging  values  using  a 
number  of  frequencies  as  low  as  M/(B/2)  =  5/10  or  0.5, 

A 

i.e.,  a  result  that  is  consistent  with  those  for  Sr  r 
for  S/N,  ■  «°  in  Fig.  6.  Although  not  shewn  in  these 
figures,  the  relatively  unaveraged  (M=l)  values  for 
show  a  peak  amplitude  of  about  1.24  at  the  frequency  for 
the  peak  ideal  value.  It  should  be  noted  that  even 

A 

though  the  values  of  at  frequencies  higher  than  Lg 
are  very  small,  their  combination  with  extremely  high  1 
values  of  the  gain  factors  at  these  higher  frequencies 

A 

yield  the  white  noise  background  to  ^  in  Figs.  7  and  8 
at  the  high  frequencies;  i.e.,  in  the  priduct  given  by  Eq. 
(94). 

A 

•  Figure  20  for  S22  and  any  S/N.^  value  shows  results  some¬ 
what  similar  to  those  for  in  the  previous  figure. 
Significant  bias  is  introduced  for  averaging  values  as 
low  as  M/(B/2)  =  5/10  or  0.5.  The  peak  value  for  S22 

at  M=1  is  about  1.22  at  the  frequency  for  the  peak  ideal 
value . 

A 

•  Figure  21  for  S  and  S/N^  =  00  indicates  results  similar 
to  those  for  sCi£i>  i.e.,  significant  bias  is  introduced 
for  averaging  values  as  low  as  M/(B/2)  =  5/10  or  0.5. 
Although  not  shown  in  the  figure,  the  relatively  unaver¬ 
aged  (M=l )  values  for  S  have  peak  values  of  about  1.6 

V  V 

at  the  L_  and  L_  frequencies  for  the  peak  ideal  values. 

h  ®2 
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A 

•  Figure  22  for  S  and  S/N.  *  10  indicates  a  significant 

yy  i 

increase  in  the  amplitude  of  Syy  relative  to  those  for 
S/t^  =  «  for  frequency  averaging  in  the  range  0.5  <  M/(B/2) 
<  2.  Although  not  shown  in  this  figure,  the  rela¬ 
tively  unaveraged  (M=l)  values  for  S  can  be  an  order 

«/  v 

of  magnitude  higher  than  the  ideal  values,  which  are 
furthermore  characterized  by  a  white  noise  structure 
across  the  whole  frequency  spectrum.  These  latter  char¬ 
acteristics  are  reflected  in  this  figure  for  averaging 
as  high  as  M/(B/2)  =  2.  The  results  of  this  figure  in- 

A 

dicate  that  the  ideal  S  would  not  be  detectable  at  any 

y  y 

amount  of  frequency  averaging  for  low  frequency  signals 
with  signal-to-instrument  noise  ratios  as  high  as  10. 

A 

•  Figure  23  for  S  and  S/N.  =  1  indicates  increases  in 

u  «y  *4' 

the  amplitude  of  the  output  auto-spectral  density  func¬ 
tion  that  are  larger  by  nearly  an  order  of  magnitude 
than  those  for  S/Ni  *  10  in  the  previous  figure.  Al¬ 
though  not  shown  in  Fig.  23,  the  relatively  unaveraged 

A 

(M=l)  values  for  S  can  now  be  about  two  orders  of 
magnitude  higher  than  the  ideal  values,  which  are  again 
characterized  by  a  white  noise  structure.  As  in  the 
previous  figure,  the  results  in  Fig.  23  Indicate  that 

A 

S  would  not  be  detectable  at  any  frequency  averaging 

y  y 

for  low  frequency  signals  with  signal-to-instrument  noise 
ratios  within  the  threshold  of  detectability. 

4.1.4  Estimated  Partial  Coherences 

Figures  24-26  show  the  estimated  partial  coherences  for  the 
gradiometer  signal  channel  (Fig.  2)  after  subtracting  out  the 
predictions  from  the  gradiometer  noise  channel  (Eq.  106).  The 
estimated  partial  coherences  are  shown  as  a  function  of  fre¬ 
quency  averaging  and  the  signal-to-instrument  noise  ratio.  The 
ideal  value  of  the  partial  coherence  for  negligible  instrument 
noise  is  unity  (Eq.  114).  This  set  of  figures  also  includes 


FIGURE  19.  Comparison  of  ideal  and  estimator  spectra  for  subsidiary  "signal 
channel  for  low-frequency  filters  with  Lg 2^ *-9 1  =  4.4. 


o 


FIGURE  22.  Comparison  of  Ideal  and  estimator  spectra  for 
gradiometer  output  for  low-frequency  filters 
with  Lg2/Lg1  *  4.4  for  S/N^  =  10. 
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FIGURE  23.  Estimator  for  gradiometer  output  spectra  for  low- 
frequency  filters  with  Lg g / Lg i  =  4.4  for  S/N^  »  1. 

Note  that  corresponding  ideal  spectra  (Figs.  21,  22) 
would  be  very  small  in  the  much  larger  scales  for 


mmBsmm 


Fig.  27  for  the  input  coherence  (Eq.  100),  which  has  an  ideal 
value  of  zero  for  uncorrelated  inputs.  These  figures  show  the 
following : 

p 

•  Figure  24  for  the  partial  coherence,  y*  2,  and  S/N1  =  ~ 
indicates  a  sharp  deviation  from  unity  near  the  L  fre- 

t=l 

quency,  a  deviation  that  is  aggravated  with  increased 
frequency  averaging  given  by  M/(B/2).  Hence,  the  par¬ 
tial  coherence  would  be  a  misleading  indicator  for  low 
frequency  signals,  since  it  shows  significant  deviations 
from  unity  even  for  M/(B/2)  =  5/10  or  0.5.  It  should  be 
noted  that  the  numerical  noise  near  Lg2  in  the  gain  fac¬ 
tor  |H1 | ,  as  indicated  in  Fig.  9  is  reflected  in  the 
partial  coherence  by  (a)  even  sharper  deviations  from 
unity  at  Lg2  than  at  Lg^  and  (b)  deviations  that  are 
independent  of  M. 

•  Figure  25  for  the  partial  coherence  and  S/Ni  =  10  exhi¬ 
bits  the  expected  impact  of  instrument  noise,  as  indi¬ 
cated  by  its  drastic  deviations  from  unity  across  the 
frequency  spectrum;  a  deviation  that  becomes  worse  with 
Increased  averaging  given  by  M/(B/2).  Although  not  shown 
in  the  figure,  the  partial  coherence  for  relatively  un¬ 
averaged  (M=l)  values  exhibits  a  white  noise  structure 
between  zero  and  unity  across  the  frequency  spectrum. 

•  Figure  26  for  the  partial  coherence  and  S/N^  *  1  indi¬ 
cates  results  similar  to  those  in  the  previous  figure 
for  S/Ni  =  10.  The  results  of  Figs.  25  and  26  suggest 
that  even  very  high  signal-to-instrument  noise  ratios 
(S/l^  >  10)  may  cause  significant  deviations  of  the 
partial  coherence  from  unity. 

2 

•  Figure  27  for  the  input  coherence  y,-  for  any  slgnal-to- 

ic  2 

noise  ratio  shows  that  the  conditions  y12  0  for  un¬ 

correlated  inputs  requires  frequency  averaging  given  by 
M/(B/2)  =  2.  Although  not  shown  in  the  figure,  the 
relatively  unaveraged  (M=l)  values  for  y^2  show  peak 
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values  near  unity,  but  with  the  white  noise  structures 
shown  in  Fig.  27.  Although  averaging  with  M/(B/2)  >  2 
might  be  desirable  to  reduce  the  input  coherence  for  the 
uncorrelated  stationary  input  series,  such  extreme  aver¬ 


aging  would  deteriorate  the  biasing  of  the  signal  S 
as  indicated  in  Figs.  7  and  8. 


Ml 


4.1.3  Overall  Results  for  Low-Frequency  Gradiometer  with 


-92— 9l- 


The  overall  results  for  the  reconstruction  of  an  imbedded 

low- frequency  signal  in  a  gradiometer  channel  characterized  by 

L„  /Le  =  4.4  indicate  the  following 
°2  °1 


For  the  ideal  condition  of  no  instrument  noise,  the  al¬ 
gorithm  introduces  significant  numerical  noise  in  the 
estimates  of  the  gradiometer  transfer  functions  at  the 
frequency  Lg2  for  the  signal  channel  and  at  Lg  for  the 
gradiometer  "noise"  channel  (Figs.  9,  14).  1 


•  For  a  signal- to- instrument  noise  ratio  of  10,  the  re¬ 
quired  number  of  frequencies  in  the  averaging  of  the 

A 

gradiometer  signal  spectrum  (S^  ^  )  and  input  coherence 
(Yj  2^  is  between  once  and  twici  ihe  bandwidth  of  the 
specified  spectra,  i.e.,  1  <  M/(B/2)  <  2  (Figs.  7,  27). 

•  For  a  signal- to- instrument  noise  ratio  of  10,  the  esti¬ 
mates  of  the  gradiometer  transfer  functions  become  un¬ 
detectable  due  to  both  the  algorithm  numerical  noise 
(Figs.  9,  14)  as  well  as  the  white  noise  at  frequencies 

higher  than  Le  (Figs.  11,  16). 

2 

•  For  a  signal- to- instrument  noise  ratio  of  10,  the  ideal 

A 

spectrum  of  the  gradiometer  output  Syy  is  undetectable 
regardless  of  the  magnitude  of  the  M/(B/2)  parameter 
(Fig.  22). 


•  The  partial  coherence  for  the  gradiometer  signal  channel 
is  not  a  useful  indicator  of  the  fidelity  of  gradiometer 
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signal  reconstruction  because  it  has  large  deviations 
from  unity  even  for  the  ideal  conditions  of  no  instrument 
noise  and  regardless  of  the  magnitude  of  the  M/(B/2)  fre¬ 
quency  averaging  parameter  (Pig.  24). 

•  For  signal- to- instrument  noise  within  the  range  of  de¬ 
tectability  (i.e.,  S/N^  =  1),  the  ideal  gradiometer  sig¬ 
nal  spectrum,  Sg  g  ,  (Fig.  2)  would  be  undetectable  be¬ 
cause  of  the  pre^ePlce  of  white  noise  comparable  to  the 

A* 

magnitude  of  the  estimated  signal  spectrum  S_  r  (Fig.  8). 

nn 

•  The  white  noise  structure  of  the  estimated  gradiometer 

signal,  Sr  c  ,  and  estimator  for  the  gradiometer  output, 

-  ^1*1 

Syy,  is  consistent  with  the  corresponding  appreciable  white 
noise  of  the  estimated  gain  factors  |H^|  and  | I  of  the 
gradiometer  transfer  functions  (Figs.  10,  15). 
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Estimated  input  coherence  for  two  subsidiary 
current  meter  channels  for  low-frequency 
filters  with  Lgg/Lg,  =  4.4. 
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4.2  LOW-FREQUENCY  GRADIOMETER  WITH  L  /L 

92  91 

Figures  28-49  show  in  the  subsequent  subsections  the  re¬ 
sults  of  the  numerical  simulations  for  Case  II  in  Fig.  5, 
which  are  characterized  by  Lg2  =  60,  Lgl  =  50,  Lg2/Lgl  =1.2 
and  B/2  =  20. 


4.2.1  Estimated  Gradiometer  Signal  Spectrum 

Figures  28-30  show  comparisons  of  ideal  and  estimated 
spectra  for  the  gradiometer  signal  channel  as  a  function  of  sig- 
nal-to-instrument  noise  ratio  for  S/N^  =  °°,  10,  1  respectively. 
These  figures  indicate  the  following 

•  Figure  28  compares  ideal  and  estimated  gradiometer  signal 
spectra  for  the  ideal  conditions  of  S/N.^  =  »  and  M/(B/2) 
values  that  are  as  high  as  the  bandwidth  of  the  specified 
spectrum.  These  results  are  similar  to  the  corresponding 
ones  in  Fig.  6  for  the  previous  case  with  Lg^/Lg^  =  4.4 
and  a  bandwidth  half  the  present  value;  i.e.,  the  esti¬ 
mated  spectrum  is  underestimated  with  about  the  same  bias 
for  M/(B/2)  =  0.5  in  both  cases.  Figure  28  shows  then  that 
an  unbiased  estimate  of  the  gradiometer  signal  spectra  re¬ 
quired  M/(B/2)  =  0.25  for  this  case. 

•  Figure  29  compares  Ideal  and  estimated  gradiometer  signal 
spectra  for  S/t^  =10.  As  in  the  corresponding  previous 
case  (Fig.  7)>  the  increased  noise  requires  a  higher 
M/(B/2)  parameter  relative  to  that  for  S/Ni  =  «  so  as 

to  obtain  an  unbiased  estimated  signal  spectrum.  However, 
the  results  in  Figs.  7  and  29  show  that  a  lower  Lg^/Lg^ 
ratio  tends  to  decrease  the  M/(B/2)  parameter  for  an 
unbiased  estimator. 


Figure  30  compares  ideal  and  estimated  gradiometer  signal 
spectra  for  S/I^  =  1.  These  results  indicate  no  relative 
effect  of  the  Lg^/Lg^  ratio,  i.e..  Fig.  8  for  Lg^/Lg^  =  4. 
and  Fig.  30  for  Lg^/Lg^  =  1.2  show  the  same  results. 


1M1-7I-3S 


FIGURE  28.  Comparison  of  Ideal  and  reconstructed  signal  spectra 
for  low-frequency  filters  with  Lgg/Lgi  =1.2  for  S/N 
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30.  Comparison  of  ideal  and  reconstructed  signal  spectra 
for  low-frequency  filters  with  Lg «/ Lg s  1.2  for  S/N 


4.2.2  Estimated  Gain  Factors 

Figures  31-35  and  36-40  show  the  estimated  gain  factors 
for  the  gradiometer  signal  and  "noise"  channels,  respectively, 
as  a  function  of  S/r^  and  M/(B/2).  These  figures  indicate  the 
following: 


Figure  31  compares  ideal  and  estimated  gain  factors  for 
the  gradiometer  signal  channel  for  S/N^  =  <*> .  The  out¬ 
standing  result  in  this  figure  is  the  disappearance  of 
the  "numerical  noise",  which  was  present  at  Lg^  for  the 


corresponding  Lg^/Lg^  =  4.4  case  (Fig. 


9).  This  result 
will  become  moreevident  from  subsequent  considerations 
of  the  respective  partial  coherences. 


•  Figure  32  shows  estimated  gain  factors  for  the  gradio¬ 
meter  signal  channel  for  S/Ni  =10.  By  comparing  the 
vertical  scales  of  Figs.  32  and  10  of  the  previous  case 
for  the  same  S/Ni  and  M/(B/2)/  it  is  apparent  that  the 
condition  L  /L  =  1.2  in  Fig.  32  has  decreased  by  a 

g  2  g 

factor  of  about  1/2  thr  white  noise  in  the  estimated 
gain  factor  for  Lg2/Lgl  =  4.4  in  Fig.  10.  These  results 
explain  the  relatively  lower  amount  of  frequency  averag¬ 
ing  required  to  estimate  Sglgl  for  S/N^  =  10  and  Lg2/ 

Lgl  =  1.2.  Figure  33  compares  ideal  and  estimated  gain 
factors  for  the  gradiometer  signal  channel  for  S/N1  =  10 
As  in  the  corresponding  cas  for  Lg2/Lgl  =  4.4  (Fig.  11) 
the  white  noise  at  higher  frequencies  in  Fig.  33  would 
prevent  detection  of  the  gain  factor  for  frequency  aver¬ 
aging  as  high  as  M/(B/2)  =  1. 

•  Figure  34  shows  estimated  gain  factors  for  the  gradiom¬ 
eter  signal  channel  for  S/N.^  =1.  By  comparing  the 
vertical  scales  of  Figs.  34  and  12  of  the  previous  case 
for  the  same  S/N.^  and  M/(B/2),  it  is  seen  that  the  con¬ 
dition  Lg2/Lgl  =  1.2  in  Fig.  34  has  decreased  by  a  fac¬ 
tor  of  about  1/3  the  noise  in  the  estimated  gain  factors 
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FIGURE  31.  Comparison  of  ideal  and  estimated  gain  factors 

for  gradiometer  signal  channel  for  low-frequency 
filters  with  Lg 2/ L g ^  =  1.2  for  S/N^  =  °°. 
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33.  Comparison  of  ideal  and  estimated  gain  factors 

for  gradiometer  signal  channel  for  low-frequency 
filters  with  Lgo/Lg,  =  1.2  for  S/N.  *  10. 
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Estimated  gain  factor  for  gradiometer  signal  channel 
for  low-frequency  filters  with  Lg0/Lg,  =  1.2  for 
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FIGURE 


35. 


Comparison  of  ideal  and  estimated  gain  factors  for 
gradiometer  noise  channel  for  low-frequency  fil¬ 
ters  with  Lg^/Lg,  =  1.2  for  S/N.  =  1. 
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for  Lg2/Lgl  =  4.1}.  Figure  35  compares  ideal  and  esti= 
mated  | H-^  |  gain  factors  for  S/N^  =  1  and  Lg^/Lg^  =  1.2. 
Figure  35  shows  that  the  ideal  |H^|  gain  factor  is  not 
detectable  due  to  the  more  than  comparable  white  noise 
background . 


Figure  36  compares  ideal  and  estimated  gain  factors  for 
the  gradiometer  "noise"  channel  and  S/l^  =  °°.  Again, 
the  outstanding  results  here  is  the  absence  of  the  numeri¬ 


cal  noise,  which  was  present  at  L 
Lg2/LSl  =  4.4  (Fig-  14 ) . 


g 


for  the  case  of 


•1 


Figure  37  shows  ideal  and  estimated  gain  factors  for  the 

gradiometer  "noise"  channel  and  S/N.^  *  10.  These  results 

for  Lg2/Lgl  ®  1.2  are  similar  to  the  corresponding  ones 

for  Lg2/Lgl  =  4.4  (Fig.  15),  except  that  the  amplitude 

of  the  white  noise  at  high  frequencies  has  increased  by 

a  factor  of  about  two  for  Lg2/Lg]=  1-2  as  compared  with 

those  for  Lg2/Lgl  *  4 . 4 ,  Figure  38  compares  ideal  and 

estimated  gain  factors  for  the  gradiometer  "noise" 

channel.  A  comparison  of  Fig.  38  with  the  corresponding 

one  for  ^g^/Lg^  =  4.4  (Fig.  16)  shows  the  absence  of 

numerical  noise  at  L~  ,  i.e.,  this  figure  also  implies 

&1 

a  persistence  of  the  numerical  noise  for  S/N^  =  00  and 

Lg2/Lgi  =  4.4  (Fig.  16)  for  S/N^  =  10  and  the  same 

L  /L  ratio. 
g2  Si 

Figure  39  shows  estimated  gain  factors  for  the  gradiometer 
"noise"  channel  and  S/N^  =1.  A  comparison  of  these  and 
corresponding  results  for  the  previous  case  (Fig.  17)  in¬ 
dicates  slmilat’  trends,  except  that  for  the  same  M/(B/2) 
value,  the  white  noise  at  high  frequency  in  Fig.  39  for 
Lg^^g^  -  1.2  is  about  twice  that  in  Fig.  17  for 
Lg,/Lgl  =  4.4.  Fig.  40  compares  ideal  and  estimated  gain 
factors  for  the  gradiomr cer  noise  channel.  As  for  S/N^  ■  10 
(Fig.  38),  the  gain  factor  would  not  be  detectable,  even 
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Comparison  of  ideal  and  estimated  gain  factors 
for  gradiometer  noise  channel  for  low-frequency 
filters  with  Lg,/Lg,  *  1.2  for  S/N<  =  ®. 


Si 


FIGURE  37.  Estimated  gain  factor  for  gradiometer  noise  chan¬ 
nel  for  low-frequency  filters  with  Lg9/Lg,  *  1.2  for 
S/Ni  =10.  6  1 
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though  the  numerical  noise  at  the  low  frequencies  in  the 
corresponding  Pig.  18  for  Lg2/Lgl  =  4.4  is  now  absent  in 
Pig.  40  for  Lg2/Lgl  =1.2. 

4.2.3  Estimator  for  Gradiometer  Output  Spectra 

Figures  41-45  show  comparisons  of  ideal  and  estimated  spec¬ 
tra  for  the  subsidiary  channels  and  gradiometer  outputs  as  a 
function  of  frequency  averaging  and  also  of  signal-to-instrument 
noise  ratio  for  the  gradiometer  output.  A  comparison  of  these 
figures  with  their  counterparts  for  the  previous  case  (Pigs. 

19-23)  shows  similar  trends  in  regard  to  the  effect  of  both 
frequency  averaging  on  the  bias  of  the  estimators  and  white 
noise  background  for  the  detectability  of  the  ideal  spectrum 
of  the  gradiometer  output  at  S/N1  other  than  infinity. 

4.2.4  Estimated  Partial  Coherences 

Figures  46-48  show  the  estimated  partial  coherences  as  a 
function  of  frequency  averaging  and  S/f^  ratio.  These  results 
are  similar  to  the  corresponding  ones  for  the  previous  case 
(Figs.  24-26),  except  that  the  larger  constant  drop  from  unity 
for  any  M  for  S/N^  =  °°  In  Fig.  24  for  Lg2/Lgl  =  4.4  is  absent 
in  the  corresponding  Fig.  46  for  Lg2/Lg-^  =  1.2  due  to  the  ab¬ 
sence  of  numerical  noise  in  the  gain  factors  for  Lg2/Lg^  =  1.2, 
(i.e..  Fig.  31  versus  Fig.  9).  The  remaining  drop  from  unity 
in  Figs.  24  and  46  for  S/N^  =  °°  also  implies  that  the  deterio¬ 
ration  of  the  partial  coherence  at  Lgl  for  Lg2/Lgl  =  4.4  is  in¬ 
dependent  of  the  numerical  noise  in  the  corresponding  transfer 
functions  (Fig.  9).  It  Is  further  noted  that  the  remaining 
deterioration  In  the  partial  coherence  remains  significant  for 
M  2  5i  i.e.,  the  averaging  required  to  obtain  nearly  unbiased 
estimates  of  the  gradiometer  signal  spectra  in  this  case  (Fig.  6). 
Finally,  Fig.  49  for  the  input  coherence  shows  trends  similar  to 
those  in  Fig.  27  for  Lg2/Lgl  =  4.4. 
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4.2.5  Overall  Results  for  Low-Frequency  Gradlometer  with  Lg^/Lg^l 


The  results  for  the  reconstruction  of  an  imbedded  low-fre¬ 
quency  signal  in  a  gradlometer  channel  characterized  by  a  ratio 
of  Lff  /L,y  ■  1.2  and  B/2  =  20,  indicate  results  similar  to  those 

&i 

for  the  previous  case  for  Ld-^/Lq.,  =  4.4  and  B/2  =  10  with  the 


following  exceptions: 


Jg2  Si 


•  For  the  ideal  conditions  of  no  instrument  noise,  the 
numerical  noise  Introduced  by  the  algorithm  in  the 
estimated  gain  factors  for  Lg^/Lg^  *  4.4  disappears 
when  Lg^/Lg^  »  1.2.  (Figs.  9  vs .  31  and  1*1  vs.  36). 

•  For  a  signal- to- instrument  noise  ratio  of  10,  the  rela¬ 
tive  window,  M/(B/2),  for  the  frequency  averaging  of 
the  estimated  gradlometer  signal  spectrum  tends  to 
decrease  from  1  <  M/(B/2)<2  for  Lg^/Lg  =  4.4  (Fig.  7) 
to  0.5  <  M/ ( B/2 ) <  1  for  Lg2/Lgi  =  1.2  Jpig.  29). 
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FIGURE  41.  Comparison  of  Ideal  and  estimator  spectra  for 
subsidiary  "signal"  channel  for  low-frequency 
filters  with  Lg,/Lg.  =  1.2. 


FIGURE  43.  Comparison  of  ideal  and  estimator  spectra  for  gradiometer  output  for 
low-frequency  filters  with  lg7/Lg,  =  1.2  for  S/N.  =  °°. 
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FIGURE  44.  Comparison  of  Ideal  and  estimator  spectra  for 
gradlometer  output  for  low-frequency  filters 
with  Lg9/Lg,  =  1.2  for  S/N,  *  10. 
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Estimated  partial  coherence  for  gradiometer 
signal  channel  for  low-frequency  filters  with 
Lg 2 / Lg i  =  1.2  for  S/N<  =  1 . 
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4.3  BROAD  BANDWIDTH  GRADIOMETER  WITH  Lg2/Lgi  = 


Figures  50-67  describe  in  the  subsequent  subsections 
results  for  case  III  in  Fig.  5,  which  is  characterized  by 
Lg2/Lgl  *  1  and  B/2  =  100. 

4.3.1  Estimated  Gradiometer  Signal  Spectrum 


the 


500 


Figures  50-52  compare  ideal  and  estimated  spectra  for  the 
gradiometer  signal  channel  as  a  function  of  frequency  averaging 
and  signal-to-instrument  noise  ratio.  For  S/N^  =  »,  Fig.  50  in¬ 
dicates  that  the  determination  of  an  unbiased  spectrum  would  re¬ 
quire  frequency  averaging  in  the  range  0.2  <  M/(B/2)  <1.  A 
comparison  of  Figs.  50-52  with  those  for  the  previous  case 
(Figs.  28-30)  indicates  a  significant  decrease  in  the  white 
noise  background  for  the  broadband  case  for  even  S/Ni  =1.  As 
a  result,  the  ideal  signal  spectrum  for  S/Ni  =  1  is  reconstructed 
almost  as  well  as  for  S/N1  =  »  for  M  =  100  or  M/(B/2)  =  1. 

4.3.2  Estimated  Gain  Factors 


Figures  53-55  and  56-58  show  the  gain  factors  for  the 
gradiometer  signal  and  "noise"  channels,  respectively.  As  ex¬ 
pected  from  the  result  in  Fig.  52,  for  the  reconstruction  of  the 
gradiometer  signal  spectrum  for  even  S/N^  =  1,  the  white  noise 
structures  at  high  frequencies  for  the  gain  factors  in  the  pre¬ 
vious  low  frequency  cases  have  now  been  reduced  drastically  even 
for  S/Ni  =1.  It  appears  then  that  the  smoothing  of  the  period- 
ogram  is  a  function  of  both  the  absolute  value  of  M  as  well  as 
the  relative  parameter  M/(B/2)  since,  as  the  bandwidth  of  the 
spectra  broadens,  the  bias  of  the  estimated  spectrum  is  reduced 
significantly  for  a  large  value  of  M. 

4.3.3  Estimators  for  Gradiometer  Output  Spectra 

Figures  59-63  show  the  reconstructed  spectra  for  the  sub¬ 
sidiary  and  gradiometer  channels  outputs.  Relative  to  the  pre¬ 
vious  low-frequency  case  with  Lg2/Lgl  =1.2  (Figs.  41-45),  the 
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outstanding  relative  result  for  this  case  is  the  adequate  recon¬ 
struction  of  the  gradiometer  output  spectra  for  S/N^  =  10  when 
M/(B/2)  =  1  (Pigs.  62  versus  44).  However,  the  reconstruction 
of  the  ideal  S  spectrum  for  S/N,  =  1  fails  to  show  a  similar 
relative  improvement  for  the  broad  bandwidth  case  (Figs.  63 
versus  45). 

4.3.4  Estimated  Partial  Coherences 

Figures  64-66  show  the  partial  coherences  for  the  broad 
bandwidth  spectra  for  signal-to-instrument  noise  ratio  S/Ni  of 
10  and  1.  For  the  required  averaging  to  detect  an  unbiased 
gradiometer  signal  spectra  for  S/N^  =  ",  0.2  <  M/(B/2)  <  1,  the 
results  in  Fig.  64  show  no  significant  difference  relative  to 
the  previous  low-frequency  case  (Figs.  24,  46),  i.e.,  the  devi¬ 
ation  of  the  partial  coherence  from  unity  for  S/Ni  =  "  remains 
essentially  independent  of  the  frequency  characteristics  of  the 
output  spectra  for  the  subsidiary  and  gradiometer  channels 
(Fig.  5).  Figure  67  for  the  input  coherence  shows  a  significant 

A 

improvement  for  broad  bandwidths,  i.e.,  y*2  +  0  as  M/(B/2)  -*•  1. 

4.3.5  Overall  Results  for  Broad  Frequency  Gradiometer  with 
L92/L9l  =  1 

The  overall  results  for  spectra  with  broad  bandwidths  indi¬ 
cate  that  the  gradiometer  signal  spectrum  and  transfer  functions 
can  be  reconstructed  adequately  for  M/(B/2)  ■  1  even  for  signal- 
to-instrument  noise  ratios  within  the  threshold  of  detectability 
(S/N^  «  1).  However,  the  spectrum  of  the  gradiometer  output 
can  be  reconstructed  only  for  signal-to-instrument  noise  ratios 
S/N.^  >  10.  The  fundamental  problem  of  the  breakdown  of  the 
partial  coherence  as  an  index  of  the  fidelity  of  reconstruction 
of  the  gradiometer  signal  spectrum  remains  Independent  of  the 
frequency  characteristics  of  such  spectra. 
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FIGURE  50.  Comparison  of  ideal  and  reconstructed  signal 
spectra  for  broad  bandwidth  filters  with 
Lgo/Lg^  =  1  for  S/N.  = 
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FIGURE  53.  Comparison  of  ideal  and  estimated  gain  factors 
for  gradiometer  signal  channel  for  broad  band¬ 
width  filters  with  Lg?/Lg,  =  1  for  S/N-  =  °°. 
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FIGURE  56.  Comparison  of  ideal  and  estimated  gain  factors 

for  gradiometer  noise  channel  for  broad  bandwidth 
filters  with  Lg 2/ Lg -j  =  1  for  S/Ni  3  00 • 
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FIGURE  57.  Comparison  of  ideal  and  estimated  gain  factors 

for  gradiometer  noise  channel  for  broad  bandwidth 
filters  with  Lg?/lg,  =  1  for  S/N,  =  10. 
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FIGURE  61 . 


Comparison  of  ideal  and  estimator  spectra  for 
gradiometer  output  for  broad  bandwidth  filters 
with  L g g / L g ^  =  1  for  S/ N H  =  «. 
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FIGURE  62.  Comparison  of  ideal  and  estimator  spectra  for 
gradiometer  output  for  broad  bandwidth  filters 
with  Lg,/Lg,  =  1  for  S/N,  -  10. 
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FIGURE  63.  Comparison  of  Ideal  and  estimator  spectra  for 
gradiometer  output  for  broad  bandwidth  filters 
with  Lg?/Lg,  =  1  for  S/N.  =  1. 
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FIGURE  64.  Estimated  partial  coherence  for  gradlometer 
signal  channel  for  broad  bandwidth  filters 
for  Lg,/Lg.  =  1  for  S/N4  =  ». 
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FIGURE  67.  Estimated  input  coherence  for  two  subsidiary 
channels  for  broad  bandwidth  filters  with 
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4.4  BROAD  BANDWIDTH  GRADIOMETER  WITH  Lg £/ Lg ]  t  1 

The  normalized  ideal  spectrum  for  the  gradiometer  output, 

S  ,  for  case  III  in  Fig.  5,  is  identical  to  those  for  Sg,g, 

u  J  J.  -i. 

and  Sg2g2  because  Lg^  =  Lg2  and  the  signal-to-gradiometer  noise 

S/Np  is  unity  (Table  1).  Since  the  results  for  this  case  with 

S/N.  =  10  for  (Fig.  51)  and  S  (Fig.  62)  indicated  that 

-*■  -L  «y  «y 

these  spectra  can  potentially  be  reconstructed  with  0.2  <  M/(B/2) 
<  1,  it  becomes  of  interest  to  investigate  if  the  algorithm 
could  do  likewise  for  other  shapes  of  the  gradiometer  spectra. 

In  order  to  test  this  possibility,  a  similar  case  was  considered 
by  displacing  Lg2  from  500  to  600  so  as  to  change  the  shape  of 
the  S  spectra  from  that  of  its  two  components  Sglgl  and  Sg2g2. 
This  modified  case  is  shown  as  case  IVa,  b,  c  in  Fig.  5  and 
Table  1,  which  included  two  subcases  so  as  to  vary  the  bandwidth 
noise  ratio,  which  were  taken  as  S/Np  =  10  and  0.1  for  the  ideal 
condition  of  no  instrument  noise  (Table  1).  These  two  values  of 
the  S/NF  ratio  would  introduce  relative  opposite  changes  in  the 
sSlSl  and  sg2S2  sPectra>  and  hence  produce  significant  changes 
in  the  amplitude  of  S  .  Figures  68-8l  describe  in  the  subse- 

*7  «/ 

quent  subsections  the  results  for  cases  IVa,  b,  c  In  Fig.  5, 
which  are  characterized  by  Lg2  =  600,  Lgl  =  500,  Lg2/Lgl  =  1.2 
and  B/2  =  100. 


4.4.1  Gradiometer  Signal -to-Noise  of  Unity 

Because  case  IVa  Is  very  similar  to  the  previous  one  with 
broad  bandwidth  and  Lg2/Lgl  =  1,  the  estimates  for  the  gradiom¬ 
eter  signal  spectra  and  transfer  functions  are  omitted.  Figures 

A 

68-70  show  the  estimator  S  for  the  case  with  a  broad  bandwidth 

«/  V 

and  Lg2/Lgl  =  1.2  as  a  function  of  frequency  averaging  and  the 
three  signal-to-instrument  noise  ratios.  A  comparison  of  the 
results  for  S/N^  =  10  (Fig.  69)  with  the  corresponding  one  for 
the  previous  case  (Fig.  62)  shows  a  better  match  in  the  shapes 

A 

of  the  ideal  (S  )  and  estimator  (S  )  spectra  for  M/(B/2)  =  1, 

«y  «y  y  %j 

as  a  result  of  the  less  pronounced  shape  of  S  caused  by  the 
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shift  of  Lg2  away  from  Lg^.  Figures  71-73  show  the  estimated 
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partial  coherences  for  the  gradiometer  signal  channel  as  a  func¬ 
tion  of  signal-to-instrument  noise  ratio,  results  that  are  simi¬ 
lar  to  those  for  the  previous  case  (Figs.  64-66). 

4.4.2  Effects  of  Gradiometer  Signal -to-Noi se  Ratio 

The  effects  of  gradiometer  signal-to-noise  ratio  are  given 
by  the  two  subcases  with  S/Np  =  10  and  S/Np  =  0.1  for  S/Ni  =  °°. 
Figures  74-77  show  results  for  S/Np  =  10  for  the  gradiometer 
signal  spectrum,  the  gradiometer  output  spectrum,  the  gain  fac¬ 
tor  for  the  gradiometer  "noise"  channel  and  the  partial  coherence 
of  the  gradiometer  signal  channel,  respectively,  whereas  Figs. 
78-81  do  likewise  for  S/Np  =  0.1. 

•  A  comparison  of  Figs.  74  and  75  for  the  gradiometer  sig¬ 
nal  and  output  spectra  for  Lg2/Lg1  =  1.2  and  S/Np  =  10 
with  those  for  Lg2/Lg1  =  1  and  S/Np  =  1  (Figs.  50,  6l) 
indicate  (a)  no  significant  difference  in  the  relative 
bias  of  the  gradiometer  signal  for  M/(B/2)  =  1,  because 
the  filter  output  Xg^(t)  in  Fig.  2  is  unaffected  Ly  the 
gradiometer  "noise"  parameter  /F,  and  (b)  no  significant 
improvement  in  the  matching  of  ideal  and  estimator  spec¬ 
tra  for  the  gradiometer  output. 

•  A  comparison  of  Fig.  76  for  the  gain  factor  of  the 
gradiometer  channel  for  S/Np  =  10  with  the  correspond¬ 
ing  one  for  the  previous  case  with  S/Np  =  1  (Fig.  56) 
shows  only  that  the  amplitude  of  the  gain  factor  for 
S/Np  =  10  has  decreased  by  a  factor  of  about  1/4  from 
that  for  S/Np  =  1. 

•  A  comparison  of  Fig.  77  for  the  estimated  partial  co¬ 
herence  for  S/Np  =  10  with  the  corresponding  one  for 

the  previous  case  with  Lg2/Lgl  =  1  and  S/Np  =  1  (Fig.  64) 
shows  no  significant  difference  in  the  deviation  of  the 
partial  coherence  from  unity  for  S/N^  =  ~. 
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•  A  comparison  of  the  results  for  S/Np  =  0.1  (Figs.  78-81) 
with  the  corresponding  ones  for  S/Np  =  10  (Figs.  74-77) 
indicate  (a)  no  significant  difference  in  the  gradiom- 
eter  signal,  (b)  a  significant  increase  in  the  amplitude 
of  the  gradiometer  spectra  by  almost  an  order  of  magni¬ 
tude  for  S/Np  =  0.1  (Fig.  79)  relative  to  that  for 
S/Np  =  10  (Fig.  75),  (c)  a  similar  corresponding  increase 
in  the  amplitude  of  the  gain  factor  for  the  gradiometer 
noise  channel  for  S/Np  =  0.1  (Fig.  80)  relative  to  that 
for  S/Np  =  10  (Fig.  76),  and  (d)  a  pronounced  enhance¬ 
ment  in  the  deviation  from  unity  in  the  estimated  par¬ 
tial  coherence  for  S/Np  =  0.1  (Fig.  8l)  as  compared  with 
that  for  S/Np  =  10  (Fig.  77). 

4.5  RECOVERY  OF  PARTIAL  COHERENCE  FROM  S/Ni  RATIO 

For  the  ideal  case  of  no  instrument  noise  (S/Ni  =  °°),  the 
fidelity  of  the  reconstruction  of  the  gradiometer  output,  y(t), 
requires  that  the  partial  coherence  be  unity  at  any  frequency 
(Eq.  114  or  46,  59-61).  However,  the  previous  results  for  S/N^  =  00 
and  M/(B/2)  ~  1  (as  required  for  detectability  of  the  gradiometer 
signal  spectra)  indicate  that  the  partial  coherence  deviates 
from  unity  independently  of  the  frequency  characteristics  of  the 
specified  spectra  (Figs.  24,  46,  64,  71)  and  for  gradiometer  sig- 
nal-to-noise  ratios  as  high  as  S/Np  =  10  (Fig.  77).  Furthermore, 
the  previous  results  for  a  signal-to-instrument  noise  ratio  that 
is  significantly  above  the  threshold  of  signal  detectability  (i.e. 
S/N^  «  10)  indicate  a  rather  drastic  collapse  of  the  partial  co¬ 
herence  away  from  unity  (Figs.  25,  47,  65,  72).  Hence,  it  becomes 
of  interest  to  determine  the  magnitude  of  the  highest  signal-to- 
instrument  noise  ratio  for  a  recovery  of  the  partial  coherence  to 
the  values  given  by  the  ideal  case  of  S/N.^  =  00 .  Figures  82-87 
show  these  results  for  the  cases  considered  in  Fig.  5.  These 
figures  indicate  the  following: 


136 


•  Figures  82  and  83  for  the  narrow-band,  low-frequency 
gradiometers  considered  in  cases  I  and  II  (Table  1), 
respectively,  show  that  signal-to-instrument  noise 
ratios  as  high  as  S/N^  =  104  fail  to  recover  the  par¬ 
tial  coherence  toward  unity. 

•  Figures  84  to  87  for  the  broad-band,  high-frequency 
gradiometers  considered  in  cases  III  and  IV  (Table  1), 
show  that  signal-to-instrument  noise  ratios  in  the  range 
102  <  S/N^  <  10 3  are  required  to  recover  the  partial  co¬ 
herence  toward  unity. 

The  foregoing  results  are  of  great  significance  for  prac¬ 
tical  applications,  since  such  extremely  high  signal-to-lnstru- 
ment  noise  ratios  would  be  required  for  the  successful  process¬ 
ing  of  time  series  data  even  in  the  case  that  it  would  become 
possible  to  extract  the  estimated  spectra  with  minimum  bias; 
i.e.,  when  the  partial  coherence  for  S/Ni  =  °°  would  indeed  be¬ 
come  unity  by  removing  the  bias  of  the  spectral  densities  in 
Eq.  (68)  for  a  multi-channel  linear  system  as  given  by  Eq.  (15) 
instead  of  the  more  general  form  of  Eq.  (14). 
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FIGURE  68.  Comparison  of  ideal  and  estimator  spectra  for 
gradiometer  output  for  broad  bandwidth  filters 
with  Lg2/Lg1  =  1.2  for  S/N.  =  ». 
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FIGURE  75.  Comparison  of  ideal  and  estimator  spectra  for 
gradiometer  output  for  broad  bandwidth  filters 
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FIGURE  77.  Estimated  partial  coherence  for  gradiometer  signal 
channel  for  broad  bandwidth  filters  with  Lg0/Lg,  * 
1.2  for  S/Nf  *  10.  £  1 
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FIGURE  81.  Estimated  partial  coherence  for  gradlometer  signal 
channel  for  broad  bandwidth  filters  with  Lg?/Lg,  » 
1.2  for  S/Nc  =  0.1.  41  1 
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FIGURE  83.  Estimated  partial  coherence  for  gradlometer  signal 
channel  for  low-frequency  filters  with  Lg„/Lg,  =  1 
for  S/N .  =  104.  2  1 
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FIGURE  84.  Estimated  partial  coherence  for  gradiometer  signal 
channel  for  broad  bandwidth  filters  with  Lg9/Lg,  ■ 
1  for  S/N 1  =  10z.  t  1 


5.  CONCLUDING  REMARKS 


A  preliminary  evaluation  of  how  well  a  mean-squared  esti¬ 
mation  algorithm  can  decompose  superposed  time  series  in  a  main 
data  channel  has  been  made  by  (a)  reviewing  the  basic  assump¬ 
tions  for  such  a  reconstruction  algorithm,  (b)  doing  the  simplest 
and  most  revealing  test  of  the  algorithm  as  applied  to  a  two- 
input  gradiometer  and  two  subsidiary  channels,  (c)  using  aibi- 
trary  but  known  spectra  to  simulate  the  gradiometer  signal  chan¬ 
nel,  so  as  to  provide  a  reference  to  test  the  corresponding 
estimated  spectra  to  be  derived  from  the  algorithm,  and  (d)  per¬ 
forming  parametric  numerical  simulations  of  the  reconstruction 
process  based  on  stationary  random  series  inputs  to  simulated 
gradiometer  and  subsidiary  channels.  These  parametric  investi¬ 
gations  allowed  study  of  the  effects  on  the  reconstruction  proc¬ 
ess  of  variables  such  as  the  frequency  characteristics  of  the 
simulated  gradiometer  and  subsidiary  channels,  the  gradiometer 
and  signal-to-instrument  noise  ratios,  and  the  number  of  fre¬ 
quencies  utilized  in  the  averaging  of  periodograms  so  as  to  re¬ 
duce  their  variance.  Some  of  the  most  important  points  brought 
out  in  these  investigations  are  as  follows: 

•  There  is  no  theoretical  basis  to  support  the  specialized 
linear  assumption  utilized  in  the  reconstruction  algo¬ 
rithm,  whereby  functional  linear  relationship  among  the 
subsidiary  channel  outputs  and  superposed  time  series 
in  the  main  channel  are  decoupled.  Hence,  the  degree 
of  confidence  in  the  fidelity  of  reconstruction  of  the 
imbedded  signal  in  the  main  data  channel  is  no  greater 
than  the  a  priori  knowledge  concerning  the  validity  of 
the  above  decoupling. 
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•  The  bias  in  the  estimators  for  the  spectral  densities 
produced  by  the  smoothing  of  their  periodograms  intro¬ 
duces  deviations  from  unity  in  the  partial  coherences 
for  ideal  conditions  of  no  instrument  noise.  This 
breakdown  of  the  partial  coherence  as  an  indicator  of 
the  fidelity  of  signal  reconstruction  is  independent 
of  the  frequency  characteristics  that  were  considered 
in  the  simulation  of  the  gradiometer  and  current  meter 
channels. 

•  The  effect  of  instrument  noise  on  the  partial  coherence 
is  a  function  of  the  frequency  characteristics  of  the 
simulated  gradiometer  and  subsidiary  data  channels.  For 
low  center  frequencies  corresponding  to  peak  amplitude 
of  the  gradiometer  signal  and  noise  spectra.  Figs.  82 
and  83  indicate  that  the  partial  coherence  remains  sig¬ 
nificantly  different  from  unity  for  signal-to-instrument 
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noise  ratio  as  high  as  10  .  However,  for  high  center 
frequencies  and  broad  bandwldths,  Figs.  84  to  87  show 
that  the  partial  coherences  recover  from  the  effect  of 
instrument  noise  at  signal-to-noise  ratios  in  the  range 
102  <  <  103. 

•  For  low  center  frequencies  corresponding  to  peak  ampli¬ 
tudes  of  the  gradiometer  signal  and  "noise"  spectra,  the 
algorithm  introduces  significant  numerical  noise  when 
the  center  frequency  of  the  gradiometer  noise  spectra  is 
significantly  larger  than  that  of  the  gradiometer  signal 
spectra.  However,  this  numerical  noise  disappears  as 
the  center  frequency  of  the  gradiometer  noise  spectra 
tends  towards  that  of  the  gradiometer  signal  spectra. 

It  is  important  to  emphasize  that  the  specified  spectra 
(Fig.  5)  used  to  test  the  algorithm  did  not  include  the  fre¬ 
quency  dependence  found  in  the  Physical  Dynamics  experiments 

_  p 

(Chapter  1),  i.e.,  Sglgl  ~  f  .  Also,  although  the  numerical 


160 


simulations  were  performed  for  a  two-input  gradiometer  and  tw 
subsidiary  channels,  the  computing  code  at  IDA  can  handle  the 
full  nine  channels  of  the  Physical  Dynamics  experiments. 
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